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I. INTRODUCTION 


The Modeling and Simulation Facility of the Laboratory for Atmos- 
pheric Sciences (GLAS) is engaged in general circulation modeling studies 
related to global atmospheric and oceanographic research. Research activities 
at the Goddard Modeling and Simulation Facility (GMSF) are organized in 
three areas: GARP/Global Weather, Climate, and Ocean/Air Interactions. 

During the period of this review, GARP/Globa! Weather research was 
directed primarily toward the development of techniques for the utilization 
and analysis of the First GARP Global Experiment (FGGE) data sets. This 
involved acquisition and data file management efforts, development of improved 
techniques for the analysis of FGGE data, observing systems research, forecast 
model development, and diagnostic studies. Observing system studies were 
concerned with the development of a GLAS TIROS N sounding retrieval system 
and preparation for the joint NOAA/NASA Advanced Moisture and Temperature 
Sounder simulation study. 

Climate research focused on the continued development of the GLAS 
general circulation model for short-range climate predictions, studies of 
the sensitivity of climate to boundary conditions, and predictability studies 
to determine which features of the global climate system can be computed de- 
terministically or statistically. Ocean/air interaction studies concentrated 
on the development of models for the prediction of upper ocean currents, 
temperatures, sea state, mixed-layer depths, and upwelling zones, and on 
studies of the interactions of the atmospheric and oceanic circulation systems 
on time scales of a month or more. 

This research review presents a synopsis of extended abstracts in 
the GARP/Global Weather, Climate and Ocean/Air Interaction studies conducted 
at the Modeling and Simulation Facility during the calendar year 1979. 

A major event during the summer 1979 was the second annual visiting 
lecture series jointly Sponsored by the Modeling and Simulation Facility and 
University of Maryland. The program included many distinguished university 
scientists and was extremely popular here at GLAS, even drawing a considerable 
attendance from the greater Washington, DC metropolitan area. Abstracts of 
the summer lectures are presented in section V. 

A list of recent publications by the staff scientists and visiting 
scientists is presented in section VI. 
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TEMPERATURE SOUNDING FROM TIROS-N 
(J. Susskind) 


TIROS-N, the current operational meteorological satellite, contains 
three passive temperature sounding instruments, HIRS2, a 20-channel infrared 
radiometer, MSU, a 4-channel microwave radiometer, and SSU, a 3-channel pres- 
sure modulated infrared radiometer sounding the upper stratosphere. NOAA/NESS 
operationally produces global atmospheric temeprature soundings by analysis of 
observations from these three instruments using statistical regression methods. 
GLAS has developed a fundamentally different approach to analysis of satellite 
temeprature sounding data, which relies more heavily on the ability to account 
for the atmospheric physics giving rise to the observations rather than the 
statistical relationships of atmospheric properties to satellite observations. 
The method is believed to be superior to that in use operationally under 
partially cloudy conditions, especially over oceans. 

The iterative scheme and cloud filtering methods used in the inver- 
sion process are closely related to those of Chahine (1970, 1974). Intrinsic 
in the use of the Chahine iterative retrieval scheme is the ability to accurately 
solve the forward problem, that is, given a guess set of atmospheric and surface 
parameters, accurately compute the corresponding observations as seen by the 
satellite instruments. The method i nvol ves iterative modification of the atmos= 
pheric parameters, starting from an initial guess, until sufficient agreement 
between computed and observed radiances is reached. 

The Radiative Transfer Equation 

Given atmospheric and surface conditions, the clear column radiances 
R-j observed by a sounding channel i can be expressed as 
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i XX s i s i x Xs 


In P 
In P B i 

C 


+ p. H. T. (P ) + 
X X XS s 


dT 


[T(P) ] 


1 d In P , <’> 


d Xn P 


where is the surface emissivity averaged over sounding channel i, B-j[T] 
is the mean Planck black-body function, averaged over channel i, of the 
temperature T, T i(P) is the mean atmospheric transmittance from pressure P 
to the top of the atmosphere and evaluated at 0, the zenith angle of the 
observation, Rj 4 is an effective atmospheric emission downward flux, pjHfX-j s 
is the reflected solar radiation in the direction of the satellite, and the 
subscript s refers to surface. The integral, taken from the surface to the 
satellite pressure P, represents the upwelling atmospheric emitted radiation, 
which is a mean value of the black-body function of atmospheric temperature 
weighted by the channel-weighting function dT^/d In P. Table 1 shows the 
channels, centers, and peak of the weighting function, and other relevant 
information, for the channels on MSU and HIRS. The current analysis does 
not employ the SSU observations. 


Table 1. HIRS2 and MSU channels frequencies and peaks of the 
weighting functions or other relevant information. 


Channel 

v(cnH ) 

Peak of dt/d In P{m) 

Peak 

of B di/d In P(inb) 

HI 

668.40 

30 



20 

H2 

679.20 

60 



50 

H3 

691.10 

100 



100 

H4 

703.60 

280 



360 

H5 

716.10 

475 



575 

H 6 

732.40 

725 



875 

H7 

748.30 

Surface 



Surface 

H 8 

897.70 

Window, 

sensitive to water vapor 

H9 

1027.90 

Window, 

sensitive to O 3 


HI 0 

1217.10 

Lower tropospheric water vapor 

Hll 

1363.70 

Middle tropospheric water 

/apor 

H12 

1484.40 

Upper tropospheric water vapor 

H13 

2190.40 

Surface 



Surface 

H14 

2212.60 

650 



Surface 

HI 5 

2240.10 

340 



675 

H16 

2276.30 

170 



425 

HI 7 

2310.70 

15 



2 

H18 

2512.00 

Window, 

sensitive to 

solar 

radiation 

H19 

2671.80 

Window, 

sensitive to 

solar 

radiation 


v (GHz) 




1 

Ml 

50.30 

Window, 

sensitive to 

surface emissivity 

M2 

53.74 

500 




M3 

54.96 

300 




M4 

57.95 

70 










The transmittance functions of the HIRS channels are taken to be 
a product of dry transmittance functions, parameterized as a function of 
temperature profile and zenith angle (Mo and Susskind, J9781, and effective 
water vapor transmittances of the formTj(P, £ " a 1 w ' K/KS ' n 1 , where W Is 
Is the column density of water vapor and a; and ni are channel -dependent 
constants (Halem and Susskind, 1977). The infrared emissivity ej is taken 
as .65 or .95 for land and water, respectively, at wavelengths greater than 
10pm, and .95 or .98 for wavelengths less than 5pm (channels between 5 and 
lOum were not used in the analysis). The effective downward flux, Rjf,is 
calculated according to Kornfield and Susskind (1977). The solar radiation 
term will be discussed later. 


The microwave transmittances are taken to be products of O2 trans- 
mlttances, calculated as a function of temperature profile and zenith angle 
according to Rosenkrantz (1975), and water vapor transmittance having the 
same form as the infrared effective water vapor transmittances. Water vapor 
retrievals using the humidity sounding channels on HIRS2 are not done at this 
time, and W, the water vapor column density is estimated at the 24-h lagged 
analysis value, as is the surface pressure P s . The microwave emissivity is 
calculated from the 50.3 GHz channels, as part of the iterative scheme, 
according to 


R - / TdT - R + T. (PJ 


It 
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where Ri Is the 50.3 GHz observed brightness temperature, T s is the iterative 
surface temperature, and T(P) is the iterative atmospheric temperature profile 
used in the calculation of the upward and downward microwave fluxes emitted 
by the atmosphere. The transmittance functions are, of course, corrected for 
temperature, water vapor, and zenith angle as described earlier, but possible 
effects of liquid water attenuation on the 50.3 GHz channel are not accounted 
for. The emissivity determined from the 50.3 GHz channel observation is used, 
together with the iterative temperature profile, to calculate brightness 
temperatures for the other MSU channels. 

Accounting for Effects of Clouds on the Infrared Observations 


The infrared radiance observed in an otherwise homogeneous field of 
view, containing partial homogeneous cloud cover 1 *, is given, by a reasonable 
approximation, as 


aR i,CLD + 


(l - a) R 


i , CLR 


(3) 


Where Ri CLD and Ri ,CLR are the radiances which would have been observed if 
the fielcl of view were completed cloudy or clear, respectively. Computation 
of clear column radiances Ri ,CLR can be done routinely as in Eq. (1), but 
computation of Ri ,cld requires accurate knowledge of the optical as well as 
meteorological properties of the cloud. It is more advantageous to be able 


4 


to account for the effects of clouds Indirectly vhan to have to model their 
radiative transfer properties. A method for doing this is employed, as pro- 
posed by Chahine (1974), using observations in adjacent fields of view using 
the assumption that both fields of view are identical (up to an accountable 
effect of zenith angle). An estimate of the clear column radiance, Ri.CLR* 
can be reconstructed from the observations according to 


VCLR ' R t,l + n ‘ R i,l ' "i.J 1 ' 0) 


where Ri i is the observation for channel i in field of view j and n is given 
by ou/fa? - «i), with ag> on. Since ais dependent only on fractional cloud cover, 
n is independent of channel and spectral region. Given n, clear column radi- 
ances can be reconstructed from the observations using Eq. (4) and the effects 
of clouds are, in principle, accounted for. 

These reconstructed clear column radiances are then used in the 
iterative temperature scheme. It is seen from Eq. (4) that large values of n. 
will tend to amplify noise in the observations and are therefore undesirable. 
As shown by Chahine (1974) and Halem et al. (1978), n can be determined as part 
of an iterative scheme according to 


P (N) - R 

n w - X Jz** , (5) 

R_ * - R„ „ 


where R 7 ^) is the computed clear column radiance for the 15wn surface channel, 
using tne N™ iterative temperature profile. The scheme will converge provided 
only 4.3pm infrared channels are used for temperature sounding in the lower 
troposphere. The rate of convergence increases with the difference between the 
surface temperature and the cloud top temperature. Under some high noise, low 
contrast conditions, divergent solutions can occur in the sense that an over- 
estimate of will cause an overestimate of the reconstructed 4.3 un clear 
column radiances which, in turn, will yield an increased lower tropospheric 
temperature, producing an increased value of R 7 ^ N + leading to an increased 
nv N + 1 ), etc. 


This situation is greatly alleviated by incorporation of a lower 
tropospheric microwave observation in the determination of n. The error in 
ni N / is a result <?f either an error in R 7 I"', due to a wrong temperature 
profile or computational uncertainties such as the effect of water vapor on 
the transmittance functions of channel 7, observational errors in R 7 ,i, or 
errors in the assumption of only one degree of nonhomogeneity in the combined 
fields of view. The error in R 7 ' N ' due to a wrong temperature profile can 
be well accounted for by adjusting the computed brightness temperature for 
channel 7 by the difference in the observed and computed microwave brightness 
temperatures for channel M2 according to 




( 6 ) 
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where Tm 2 and / N \ t M 2 are the observed and calculated microwave brightness tem- 
peratures, T 7 ' N ' is the calculated clear column brightness temperature for 
channel 7, and T7 is the corrected clear column brightness temperature for 
channel 7. The corrected clear column radiance for channel 7, to be used in 
Eq. (5), is then given by 



i ,r T <"> 
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This procedure not only speeds up convergence under all conditions but also 
stabilizes the solution in the sense that an increase in the iterative tem- 
perature profile in the lower troposphere will not, to a first approximation, 
cause an increase in rp 


Determination of Surface Temperature 

Given the clear column radiances for the three window chan- 

nels 8, 18, and 19 are reconstructed according to Eq. (4). All three channels 
are relative atmospheric windows and are sensitive primarily to be surface 
(ground) temperature. The two 3.7pm channels have the advantage of being more 
sensitive to surface temperature and less sensitive to, uncertainties in surface 
emissivity and atmospheric water vapor than the 11pm window channel. They 
have the disadvantage of being affected by solar radiation during the day, 
which must be accounted for before accurate surface temperatures can be calcu- 
lated. At night, surface temperatures are taken to be the average of the 
surface temperature as determined from channels 18 and 19, where 



S, X 1 
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In general, T s j8 ar!£l Tgjg agree with each other to 1° but differ by a 
larger amount from the surface temperature as determined from the 11 pm 
window channel 8, especially when the water vapor column density along the 
path of observation is greater than 3 gm/cm 2 . 

During the day, the effects of solar radiation on the 3.7pm channels 
must be accounted for in obtaining accurate surface temperature retrievals 
from these channels. This can be done directly by subtracting PjHj x 1S (P S ) 
from and substituting the result into Eq. (8). Hj'r- (P $ ), the mean solar 

radiation across the channel, transversing the path from the sun to the earth 
and back to the satellite, can be well estimated as 2.16rr x 10" 5 [5600K]cos 

0 H T i(Ps,^H + ®) where is the solar zenith angle and the transmittance 

is computed at an effective zenith angle given by the sum of the solar and the 
satellite zenith angles. 

The danger in such a procedure is the uncertainty in Pj. If the 
surface is Lambertian and the emissivity is known, Pj , the directional reflec- 
tance, is equal to (1 - £j)/Tr. Significant errors of up to a factor of 2 can be 


t 


-~T“ 


6 


made in the above estimation of pj, which may produce errors of up to 10K in 
retrieved surface temperature. These errors arise from uncertainties inej and 
the non-Lambertian character of the surface. However, the same uncertainties 
in ej do not appreciably effect the calculated thermal radiation. 


An error in Pi will produce different errors in the surface tempera- 
ture as retrieved from the two 3.7ym channels because dT/dB in channel 19 is 
twice as great as in channel 18. Consequently, agreement to within .(L 0 . of 
T c it* and T. iq , obtained by subtracting (1 - e^/irHf Ti ( p s ) from in 

Eq. ?8), istaken as evidence of an accurate estimate value of Pi* and the 
average of T s 18 and T s 19 is taken as the surface temperature. If this 
agreement is not obtained, T s 8* obtained from the llym window channel, is 
taken as the ground temperature, provided the water vapor column density along 


vriv. ~ » « rt 

the observing path is less than 3 gm/cnr. 


As a final alternative, an attempt to solve simultaneously for T s 
and Pj is made by linearizing Eq. (1) about T s ,8> which is expected to be a 
reasonable approximation of T s , according to 


r. = eb. [t 
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s,8 
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T. (PJ + £ 

1 S 
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s,8 


(T 


T s,8> 
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Eq. (9) gives two linear equations, one for channel 18 and one for channel 19, 
for the two unknowns, T s and p, and can be solved in a straightforward irr' ier. 


Iterative Procedure 


Given an initial or guess temperature profile T^(P), and surface 
temperature T s , clear column radiances can be calculated from Eqs. (1) and (2) 
and algorithms to calculate transmittance functions and downward fuxes., .Eqs. 
(4), (5), and (7) are then used to reconstruct clear column radiances for 
the infrared channels. T $ N+T is obtained from Eqs. (8) or (9). T 1 ’* ( P ) is 
obtained in a relaxation manner analogous to that of Chahine. Chahine assigns 
a pressure Pi to each temperature sounding channel i, given approximately by 
the peak of the weighting function for channel i, representing that portion 
of the atmosphere in which local changes of temperature will have the largest 
effect on the observed radiance for that channel. Chahine uses the iterative 
equation, 


r N+l . . i 

B i [T (P^l 

B i [*r N ( p i> ] 


; (n) 
i 



(10) 


to determine a new estimate of temperature at each of the pressures P^. 
Temperature at other pressures can be obtained by a variety of interpolation 
schemes or other constraints on the solution. Given T N+I (P) and T^ * the 
entire iterative procedure can then be repeated. 
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A more convenient form of the iterative equation for temperature can 
be written in terms of clear column brightness temperatures for each 
channel, that is, the temperature of a black body with radiance R-j. The 
analogous equation to Eq. (10) becomes 


t n+1 (p ;L ) * t n (p ;L ) + t! n) 


_ T W 
i 


(ID 


Eq. (11) reflects the very good approximation that, given a temperature profile 
and brightness temperature, if the entire profile and surface temperatures 
were changed slightly by AT, the brightness temperature will also change by 
AT regardless of the frequency region or nature of the weighting function. 

The iterative scheme based on either Eq. (10) or, equivalently, 
Eq. (11) sometimes produces divergent solutions if noise in adjacent channels 
is large and of opposite sign. It was found that stability was increased, 
with only slight loss in speed of convergence, if the temperature at pressure 
Pj was modified according to a weighted difference of observed and calculated 
brightness temperatures for all channels according to 
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where Wjj is the relative change in brightness temperature for channel j 
produced by a change in temperature at pressure P^. 

From Eq. (1), to a good approximation, 
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Wij, the relative change in brightness temperature for channel j produced 
by a change in temperature at P-j , is then given by 
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Using appropriate approximate forms for the black-body function at 
microwave and infrared wavelengths, one gets 


w. . = 
13 
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for infrared channels and 
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for microwave channels. 

Fig. 1 shows the weighting ^unctions as defined in Eq. (15) for the 
six TIROS-N channels used in the analysis of the data, for U.S. standard atmos- 
phere under nadir viewing. The observations are most sensitive to atmospheric 
temperature changes below 500 mb and between 150 mb and 50 mb, with relative 
weakness in between. The six characteristic pressures used in the analysis 
are also shown in the figure. 

Gridding of Satellite Data 


In order to facilitate analysis of global data on a twice-daily basis, 
the HIRS data are gridded (averaged) prior to analysis so that one sounding is 
performed every 250 x 250 km. The 250 x 250 km area is first divided into 
four equal quadrants. In each quadrant, the HIRS spots (20 km x 20km at nadir) 
are divided into equal groups, representing the two fields of view, according 
to the brightness temperature of channel 8, the 11pm window, with field of 
view 1 containing those spots with the largest radiances for this channel. 
The radiances for all spots in a field are averaged together for each channel 
to give the radiances for each field of view to be used in the analysis of the 
data. Each field of view is assigned a zenith angle weighted according to the 
cosine of the zenith angles of the individual spots. The microwave brightness 
temperatures and zenith angle of the MSU spot (110km x 110 km at nadir) most 
overlapping the quadrant, are assigned to each field of view. Using this 
gridding scheme, approximately 17000 retrievals are performed daily. The 
computing time for this is roughly 1 h CPU time on the Amdahl V/460. 

Results 

Global retrievals were run for the period 4 January-5 February 1979 
using a 24-h lagged analysis initial guess. A shape-preserving interpolation 
scheme (Chahine, 1970), in which t( n+ ')(P) - tN ( P) was taken to be linear in 
log P, was used in analysis of the data. In all of the cases, a sounding was 
rejected if, after five iterations, the RMS difference of computed and recon- 
structed brightness temperatures for the six temperature sounding channels was 
greater than 1° or the computed and observed brightness temperature for chan- 
nel M2 differed by more than 1.5°. Approximately 30 percent of all soundings 
were rejected by these criteria, and of these, half which were rejected were 
considered clear. This indicates that only in approximately 15 percent of the 
cases was it too cloudy to perform retrievals using the infrared tropospheric 
channels. In contrast, the operational retrievals produced by NESS perform 
infrared tropospheric retrievals only 50 percent of the time. In the remainder 
of the cases, retrievals are performed using the HIRS stratospheric channels 
together with the MSU channels. These have been shown to be considerably less 
accurate compared with radiosondes. 
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Fig. 2 shows the RMS errors of layer mean temperatures from both 
the QLAS and operational retrievals, as compared to oceanic radiosondes 
colocated to within 110 km in space and 6 h in time. In the first graph, 
errors of all retrievals from each set, colocated with radiosondes are compared 
with each other. In the second graph, comparison is made only between those 
retrievals performed by NOAA under partially cloudy conditions and successful 
GLAS retrievals colocated to the same radiosondes. The GLAS retrievals are 
clearly significantly better in the troposphere under these conditions, in 
which the combination of the cloud problem and poor statistical representation 
of oceanic radiosondes in the regression data base used by NESS significantly 
degrade the operational results. 

Fig. 3 shows a comparison of the 1000-500 mb thickness fields off 
western Europe at 0300 GMT 4 January 1979, derived on the McIDAS from the GLAS 
soundings and the NESS enhanced soundings (see Atlas et al . , 1980). The 
general agreement between the fields is good, especially in the region over 
western Europe validated by radiosonde reports. The differences and their 
sources are being evaluated, particularly in regions where the enhanced 
soundings were microwave only, such as over England. 
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OCEAN SOUNDING ERROR STRUCTURE JAN 3-FEB 5, 1979 


- NESS ALL SOUNDINGS 
2.23° RMS ERROR 

- GLAS ALL SOUNDINGS 
2.36° RMS ERROR 


NESS N* OR MSU 

SOUNDINGS 
2.67° RMS ERROR 

GLAS COLOCATED 

SOUNDINGS 
2.41° RMS ERROR 


.Til 


RMS ERROR °K 

RMS errors of GLAS and NESS operational soundings over oceans 
as compared to radiosondes colocated to 110 km. a) shows all 
NESS soundings colocated to radiosondes and all GLAS soundings 
colocated to radiosondes; b) shows statistics for all NESS 
retrievals under partially cloudy conditions colocated to both 
radiosondes and GLAS retrievals as well as statistics for the 
colocated GLAS retrievals. 
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1000-500 mb thickness analyses based on enhanced TIROS-N re- 
trievals (solid lines, see Atlas et al . , 1980) and 6LAS 

TIROS-N retrievals (dashed lines). The 6 LAS retrievals were 
unedited while the enhanced retrievals were subjectively edited 
by man-computer interaction. 


ENHANCEMENT OF TIROS-N SOUNDING AND CLOUD MOTION 
WIND DATA FOR THE FGGE SPECIAL EFFORT 

(R. Atlas, G. Cole, R. Rosenberg, S. Palm, and A. Pursch) 


A "Special Effort" to produce high quality edited and enhanced 
global data sets is being conducted for the two special observing periods of 
the First GARP Global Experiment (FGGE), Greaves et al. (1979). The Special 
Effort is a joint project between NASA, NOAA, and the' University of Wisconsin. 
Man-computer Interactive Data Access System (McIDAS) terminals, developed by 
the Space Science Engineering Center (SSEC) of the University of Wisconsin, 
are being utilized by experienced meteorologists at the National Meteoro- 
logical Center (NMC) for data evaluation and quality assessment, and at the 
Goddard Laboratory for Atmospheric Sciences (GLAS) for satellite data enhance- 
ment. Interactive procedures for sounding and cloud motion wind enhancement 
were originally developed by the National Environmental Satellite Service and 
SSEC, and are described by Smith et al . (1978) and Menzel et al . (1978). More 
than 60 cases have been selected for enhancement, on the basis of meteoro- 
logical interest, by an ad hoc international committee of participating 
scientists. These cases include situations of blocking, cutoff low develop- 
ment, cyclogenesis, and tropical circulations. 

The sounding data enhancement is aimed at supplementing the opera- 
tional satellite sounding data set with higher resolution soundings in 
meteorologically active regions, and with new soundings where data voids or 
soundings of questionable quality exist. The cloud motion wind enhancement 
is aimed at the rederivation of cloud heights in suspect areas using radiance 
information from polar orbiting satellites. After the enhancement process 
is completed, a final editing and quality assessment of the enhanced data is 
performed. The data will be archived as Level 1 1 b and will be utilized in 
Level Illb analyses. 

The sounding enhancement process at GLAS is performed as fol»yWs: 
Once an area has been selected, the McIDAS operator begins by displaying all 
available conventional and special FGGE data for that area, the Level III 
analysis, and visible, infrared, and microwave images from Television Infrared 
Observational Satellite (TIROS) -N . At this stage, the operator notes where 
data deficiencies exist and, from the TIROS-N images, determines the extent 
of cloudiness and where the most intense atmospheric thermal gradients are 
located. High resolution temperature retrievals are then generated for the 
area. This is followed by a comprehensive manual editing of the retrievals 
to remove small-scale discontinuities due to cloud-induced noise while 
retaining significant meteorological structures. 

Sounding data enhancement operations began on 3 December 1979. The 
first case selected for enhancement was a strong cyclone centered southeast 
of England at 0400 GMT 4 January 1979. The 1000 mb analysis for the enhanced 
region, from conventional data at 0000 GMT 4 January 1979, is shown in Fig. 1. 
Fig. , resents the 1000-500 mb thickness analysis utilizing the 29 radiosonde 
reports <at were available for this region (with no satellite soundings) as 
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solid lines and the 1000-500 mb thickness analysis utilizing 141 enhanced 
satellite soundings (with no radiosonde data) as dashed lines. Good agreement 
between the pure SAT and pure NOSAT thickness analyses can be seen over western 
Europe, while poorer agreement is evident over England where only microwave 
retrievals could be generated, and over oceanic areas where virtually no radio- 
sonde data was available and the NOSAT analysis is questionable. 
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THE COMBINED INFLUENCE OF SATELLITE TEMPERATURE SOUNDING DATA AND 
INCREASED HORIZONTAL RESOLUTION ON GLAS MODEL FORECASTS 

(R. Atlas, M. Halern, and M. Ghi 1 ) 


Since the introduction of satellite-derived temperature soundings into 
the operational data base in 1972, a number of experiments have been conducted 
(Atkins and Jones, 1975; Desmarais et al. , 1978; Druyan et al., 1978, Ghil 
et al. , 1979a, b) to assess their influence on numerical weather forecasting* 
Object i ve measures of forecast accuracy have served as the primary tool for 
evaluating the predictive impact of the satellite temperature sounding data, 
along with some subjective evaluations by weather forecasters. The objective 
measures have generally shown either negligible or small beneficial impacts 
when averaged over many cases and continental regions. Even when specific 
forecasts show substantial objective impacts, subsequent subjective evaluations 
have shown the improvements to be of marginal value for forecasters. A principal 
reason may have been the fact that these statistical improvements occur in the 
48-to-72 h forecast range when the model forecast accuracies have themselves 
deteriorated. 

Thus, the question arises whether satellite sounding data might have 
greater influence, beneficial or otherwise, if numerical weather prediction 
(NWP) models had greater forecast skill. Recent studies by Robert (1976), 
Williamson (1978), and Quirk and Atlas (1977) have shown that forecast accuracy 
is generally improved with increased model resolution. The study here will 
report on comparisons of temperature sounding data impacts for the GLAS model 
with two different grid resolutions: a coarse 4° lat. x 5° long, grid and a 

finer 2.5° lat. x 3° long. grid. 

The influence of satellite temperature sounding data and increased 
horizontal resolution on GLAS model forecasts was evaluated by means of objec- 
tive measures of forecast accuracy and subjective comparisons of prognostic 
charts for eleven different forecasts from the February 1976 Data Systems Test 
(DST) period. Table 1 (adapted from Ghil et al . , 1979b) summarizes the effects 
of satellite data and model resolution for the same quantities and regions of 
verification, previously presented by Ghil et al. (1979a). 

Table 1. Objective assessments of forecast impact. 


Experiment 

Control 

Impact 

Si 

(%) 

RMS 

Statistical significance 
(average/standard error) 
Sj RMS 

SAT coarse 

NOSAT coarse 

4.46 

12.41 

1.85 

2.59 

NOSAT fine 

NOSAT coarse 

5.86 

7.30 

1.62 

1.35 

SAT fine 

NOSAT fine 

4.49 

7.98 

1.40 

1.71 

SAT fine 

NOSAT coarse 

10.01 

14.47 

2.67 

2.47 

L 


18 





In the table, the results of the SAT experiment with the fine model were com- 
pared with those of the coarse NOSAT experiment and with those of the fine 
NOSAT experiment. The fine NOSAT experiment appears as a control and in a 
comparison versus the coarse NOSAT. It can be seen that the effects of 
increased resolution and satellite data both act in the direction of improving 
the forecast accuracy. The two effects appear to be nearly additive, as 
indicated by a comparison of the results in the last three rows of the table. 
In addition, the SAT experiment, with the fine model, consistently gave impacts 
that are significantly higher than those of the corresponding experiment with 
the coarse model, when measured against the common standard of the coarse 
NOSAT experiment (viz., the second and the last row of the table). 

Tables 2 and 3 present the consensus of three experienced forecasters 
in their subjective assessment of the prognostic charts. During the evaluation, 
the forecasters had no knowledge of which prognostic charts were made from SAT 
or NOSAT initial conditions. Following Atkins and Jones (1975), each forecast 
was classified according to the following scale: 

A - SAT signficantly better than NOSAT, 

B - SAT better than NOSAT, 

C - SAT and NOSAT of equal quality, 

D - SAT worse than NOSAT, 

E - SAT significantly worse than NOSAT. 

Each forecaster performed his evaluation independently. This was followed by 
a discussion of the prognostic differences, after which the ratings were 
averaged to arrive at the consensus. 

The results of the comparisons for the sea level pressure and 500 mb 
height fields generated from the GLAS 4° 1 at . x 5° long, model and verified 
over North America are presented in Table 1 (reproduced from Ghil et al . , 
1979a), From this table it can be seen that at sea level none of the prog- 
nostic charts have been classified a< A or E and that an almost equal number 

of B and D marks have been assigned. Although no significant impact has 
occurred at this atmospheric level, a tendency for more beneficial impacts 
at 48, 60, and 72 h is evident. At 500 mb, two 72-h prognostic charts have 

been classified as A and there were no E marks. In addition, more than three 

times as many B than 0 marks have been assigned. Impacts at this level are 
most consistently beneficial at 60 and 72 h. 

Table 2 presents the results of the synoptic evaluation of sounding 
data impact using the GLAS 2.5° x 3° model. This table shows that a more 
substantial beneficial impact from satellite sounding data has resulted from 
the use of the higher resolution model. No significant negative impacts have 
occurred while five sea level pressure and four 500 mb prognostic charts have 
been significantly improved. Once again the effect of satellite data is 
largest during the latter half of the forecast period. 
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Table 2. Results of a subjective assessment of SAT and NOSAT forecasts 
using the GLAS 4° x 5° model. 


Category 

12 h 

24 h 

36 h 

48 h 

60 h 

72 h 

! 

Total 




Sea level 

pressure 




. 

A 

0 

0 

0 

0 

0 

0 

0 

B 

0 

1 

1 

4 

3 

4 

13 

C 

8 

8 

7 

5 

6 

5 

39 

D 

3 

2 

3 

2 

2 

2 

14 

E 

0 

0 

0 

0 

0 

0 

0 



500 mb geopotential height 




A 

0 

0 

0 

0 

0 

2 

2 

B 

2 

2 

2 

3 

4 

4 


C 

8 

9 

9 

5 

7 

4 


D 

1 

0 

0 

3 

0 

1 


E 

0 

0 

0 

0 

0 

0 

0 ! 












Two cases of significant forecast impact, the forecasts from 19 and 
11 February 1976, were selected for a more detailed subjective analysis. The 
72-h sea level pressure prognoses from these cases are presented here as illus- 
trations of the combined influence of satellite data and increased resolution. 

72-h Forecast from 0000 GMT 19 February 1976 

This was a case in which a significant improvement to the predicted 
displacement of an intense winter storm occurred in the latter part of the 
forecast period. The initial conditions for this forecast showed a moderately 
intense low pressure system, associated with an upper level short wave trough, 
located off the northwest coast of the U.S. As the system moved inland, a new 
low developed along an already existing stationary front and became the dominant 
feature by 1200 GMT on 19 February. During the next 24 h, the low moved south- 
eastward and intensified, after which time it recurved and then accelerated 
toward the northeast. This cyclone was accompanied by heavy snow, blizzard, or 
near-blizzard conditions in Colorado, Kansas, Nebraska, Iowa, Michigan, and 
Wisconsin, and tornadoes or severe thunderstorms in Kansas, Oklahoma, Texas, 
Arkansas, Missouri , Illinois, Louisiana, and Mississippi. Most of the severe 
weather occurred toward the end of the forecast period after the recurvature of 
the cyclone to the northeast. 

Figs, la and lb depict the 72-h SAT and NOSAT sea level pressure 
prognoses for the 4° x 5° model, while the corresponding SAT and NOSAT prognoses 
for the 2.5° x 3° model are shown in Figs. 1c and Id. The verifying analysis 
is shown in Fig. le. Comparison of these charts reveals little difference 
between the two 4° x 5° model forecasts, however, a very significant improvement 
in the predicted displacement of the surface low has resulted from the inclusion 
of satellite data in the 2.5° x 3° model. The high resolution SAT prediction 
indicated recurvature of the low to the northeast while the other predictions 
did not. As a result, there was a more than 1000 km reduction in the vector 
error of the predicted position. 

Since most of the severe weather occurred in conjunction with the 
recurvature of the low, it is clear that significantly improved local weather 
forecasts could have resulted from the use of the high resolution SAT prognoses 
in this situation. To investigate this point further, the high resolution SAT 
and NOSAT predictions of local convective instability, and also t'ne destabiliza- 
tion by differential equivalent potential temperature advection, were compared. 
It was found that there was a 51 percent reduction in the RMS errors of con- 
vective instability and a 22 percent reduction in RMS errors of differential 
advection, in the prediction which included satellite-sounding data. 

72-h Forecast from 0000 GMT 11 February 1976 

In this case, a weak cyclone formed along a stationary front in 
southwest Canada, moved southeastward while intensifying during the first 40 h 
of the period, and then recurved to the east northeast. 

Figs. 2a and 2b depict the 72-h SAT and NOSAT sea level pressure 
prognoses for the 4° x 5° model, while the corresponding SAT and NOSAT prognoses 
for the 2.5° x 3° model are presented in Figs. 2c and 2d. Fig. 2e shows the 
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j. 1. The sea level pressure fields on 0000 GMT 22 February 1976: 

a. the 72-h forecast from the 4° x 5° GLAS model with SAT 
initial conditions, 

b. the corresponding NOSAT forecast using the 4° x 5° model, 

c. the corresponding SAT forecast using the 2.5° x 3° GLAS 
model ; 

d. the corresponding NOSAT forecast using the 2.5° x 3° GLAS 
model , 

e. the verifying NMC analysis. Dots represent past positions 
of the cyclone center at 12-h intervals beginning 0000 GMT 
19 February 1976. 
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The sea level pressure fields on 0000 GMT 14 February 1976: 

a. the 72-h forecast from the 4° x 5° model with SAT initial 
conditions, 

b. the corresponding NOSAT forecast using the 4° x 5° GLAS 
model ; 

c. the corresponding SAT forecast using the 2.5° x 3° GLAS 
model ; 

d. the corresponding NOSAT forecast using the 2.5° x 3° GLAS 
model ; 

e. the verifying NMC analysis. 
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verifying analysis. Comparison of these figures reveals that an improvement in 
the predicted displacement of the surface low has resulted from the inclusion 
of satellite data in the 4° x 5° model, but a much larger improvement has 
occurred in the 2.5° x 3° model. The high resolution SAT prognosis (Fig. 2c) 
is substantially better than the NOSAT (Fig. 2d) in its prediction of the 
cyclonic circulation over southeastern Canada and the northeastern U.S., as 
well as the anticyclone near the Great Lakes and the pressure trough which 
extends southeastward from southwest Canada. 
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A CASE STUDY OF THE GROWTH OF PROGNOSTIC DIFFERENCES BETWEEN 
GLAS MODEL FORECASTS FROM SAT AND NOSAT INITIAL CONDITIONS 

(R. Atlas, R. Rosenberg, and S. Palm) 


The high .resolution GLAS model forecasts from OOOO GMT 19 February 
1976, presented by Atlas et al . (1980), represent a case in which a signifi- 
cant beneficial impact resulted from the assimilation of satellite sounding 
data. Comparison of the 72-h forecasts from SAT and NOSAT initial conditions 
(Figs, lc and Id of Atlas et al . , 1980) revealed that the predicted displace- 
ment of an intense surface low was similar for the first 36 to 48 h of the 
forecast, but then diverged rapidly. A study has been conducted to determine: 
(1) why the surface lows in the two forecasts suddenly diverge and move 
differently during the latter half of the forecast and (2) if differences in 
the upper level forcing for the cyclone at the time the two forecasts diverge 
can be traced to specific initial state differences between the SAT and NOSAT 
systems. 


The initial sea level pressure and 1000-500 mb thickness, and 
300 mb height and vorticity analyses, for the SAT and NOSAT forecasts from 
19 February are displayed in Figs. 1 and 2. Comparison of these figures 
reveals no substantial difference in the initial representation of the surface 
low located at 46N, 132W, off the northwest coast of the U.S., although a 
modification of the 1000-500 mb thickness pattern has occurred. This results 
in an enhanced variation of thermal vorticity and thermal advection across the 
low center in the SAT case. At 300 mb, the assimilation of satellite sounding 
data has resulted in an intensification of the upper level vorticity maximum 
associated with the surface low and an elimination of the vorticity trough 
extending southward from this maximum. The satellite data has increased the 
northerly component of the wind and the cyclonic shear of the westerly com- 
ponent of the wind from 40N to 47N to the immediate west of the upper level 
trough. It also increased the anticyclonic shear of the westerly wind below 
40N (wind component charts not shown). This results in a northward shift of 
the upper level vorticity advection areas in this region, such that there is 
stronger positive vorticity advection and positive vorticity tendencies to 
the east of the 300 mb vorticity maximum, and stronger negative vorticity 
advection and negative vorticity tendencies to the west of the maximum in the 
SAT case. This modification of the vorticity advection by satellite sounding 
data was found to be in agreement with satellite cloud imagery (see Hales, 
1979, for discussion of technique). Additional modifications to the initial 
state by the assimilation of satellite data have been considered in this study 
but will not be described here. 

Figs. 3-10 display the evolution of the sea level pressure, 1000- 
500 mb thickness, and 300 mb height and vorticity patterns at 12-h intervals 
for the first 48 h of both the SAT and NOSAT predictions. During the first 
12 h of the forecast, both SAT and NCSAT move the surface low inland to 46N, 
123W (Fig. 3). The SAT surface low has intensified while the NOSAT low has 
filled slightly. However, the SAT thickness gradient has weakened relative 
to the NOSAT, and there is now a larger variation of thermal vorticity and 
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stronger thermal vorticity advection across the low in the NOSAT case. At 
300 mb, a small difference in the movement of the height trough and vorticity 
maximum has occurred (Fig. 4). The NOSAT vorticity maximum is located at 
44N, 128W, while the SAT vorticity maximum has moved about 2.5° longitude 
further downstream. A comparison of the vorticity advection patterns and 
tendencies at 2-h intervals throughout this period indicated that the initial 
tendency for greater movement in the SAT system has been maintained. 

During the next 24 h of the forecast, this trend continues and 
slightly larger differences between the SAT and NOSAT systems begin to develop. 
At 0000 GMT 20 February (Figs. 5 and 6), the NOSAT system has forecast the 
surface low to deepen slightly and move to 41. 3N, 108W, and the 300 mb vor- 
ticity maximum to move to 42N, 123W. SAT forecasts a less organized surface 
low with three centers evident at this time but has maintained a somewhat 
stronger pressure gradient than the NOSAT to the west and southwest of the 
southernmost low center. SAT forecasts the 300 mb vorticity maximum to move 
to 41 N, U7W, 5° longitude further downstream than the NOSAT maximum. 
However, in both forecasts, the strongest positive vorticity advection is to 
the southwest of the surface low position. By 1200 GMT 20 February, the 
NOSAT surface low is located at 41. 2N, 105W, while SAT forecasts two low 
centers: one at 41. 2N, 102W; the other at 36. 2N, 102W (Fig. 7). Substantial 

differences in the low level thermal advection patterns are also evident at 
this time. In the NOSAT case, weak cold advection exists directly behind the 
surface low, with the strongest cold advection located further southwest 
beneath and slightly to the west of the upper level trough. This contributes 
to the deepening and slow rate of movement of the upper level trough in the 
NOSAT prediction. SAT predicts substantially stronger cold advection directly 
behind the surface low. 

The relationship of low level thermal advection to the movement of 
surface lows is well known (Sutcliffe, 1947; Petterssen, 1954). Warm advec- 
tion contributes to falling pressure and the generation of cyclonic vorticity 
in advance of cyclones while cold advection contributes to rising pressure 
and the destruction of cyclonic vorticity to the rear of cyclones. Thus, 
the gradient of thickness advection across a cyclone center is an important 
indicator of the rate of movement of that cyclone. In this case, the dif- 
ferences in the gradient of thickness advection that are established after 
36 h are primarily associated with the differing rates of movement of the 
upper level vorticity maximum. Because of the different movements of the 
vorticity maximum relative to the surface low, there was an intensification 
of the pressure gradient to the southwest of the low center in the SAT 
forecast, while a broad area of surface pressure falls associated with the 
positive vorticity advection aloft weakened the pressure gradient to the 
immediate southwest of the surface low in the NOSAT forecast. 

This effect amplified during the next 12 h in such a way that by 
0000 GMT 21 February the gradient of thickness advection across the low center 
is nearly twice as strong in the SAT forecast as in the NOSAT. At this time, 
the surface low is located at 36. 5N, 102W, in the NOSAT and 36N, 99W, in the 
SAT (Fig. 9). The 300 mb vorticity maximum is located at 36. 2N, 111W, in the 
NOSAT and 38N, 102W, in the SAT (Fig. 10). The differing phase relationships 
between the upper level vorticity maximum and the surface low coupled with 


27 















Fig. 9. As in Fig. 1 except for 0000 GMT Fig. 10. As in Fig. 2 except for 0000 GMT 
21 February 1976. 21 February 1976. 


the differing thickness advection patterns result in the diverging paths of 
the surface low throughout the remainder of the forecast. In the SAT case, 
there is strong positive vorticity advection and warm advection to the east 
and northeast of the surface low and weaker vorticity advection coupled with 
strong cold advection to the immediate west and southwest of the surface low. 
This results in a strong isallobaric gradient across the low center such that 
recurvature to the east northeast occurs. In the NOSAT case, weak positive 
vorticity advection is coupled with warm advection to the northeast of the 
surface low while strong positive vorticity advection is coupled with weak 
cold advection to the immediate west and southwest of the surface low. As a 
result, the NOSAT surface low becomes "locked in" and does not progress 
eastward after this time. 
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SYNOPTIC EVALUATION OF GLAS AND NMC HIGH RESOLUTION 
FORECASTS FROM 19 AND 11 FEBRUARY 1976 

{R. Atlas) 


The synoptic evaluation of the impact of satellite temperature 
sounding data on eleven high (-280 km) resolution GLAS model forecasts 
(Atlas, et al. , 1980) revealed no examples of significant negative impacts 
and two cases in which major improvements to sea level pressure forecasts 
occurred: the 72-h forecasts from 0000 GMT 19 and 11 February 1976. In 
this paper, the GLAS SAT predictions for both of these cases are compared 
with the corresponding NOSAT forecasts from the National Meteorological 
Center (NMC) 7-level model. This model, which has a horizontal resolution 
of 190.5 km, has shown considerable improvement over the coarser mesh 6-level 
model and has been adopted by NMC as their operational forecasting model. 
Comparison of the GLAS SAT forecasts with the NMC high resolution NOSAT 
forecasts should provide additional evidence concerning the meteorological 
significance of the GLAS sounding data impact results. 

Forecasts from 0000 GMT 19 February 1976 

Figs. 1-6 depict the NMC 7-level NOSAT and GLAS SAT 72-h sea level 
pressure (solid lines) and 1000-500 mb thickness (dashed lines) forecasts 
and corresponding verification at 12-h intervals. Comparison of the two 
forecasts reveals that the GLAS SAT prediction of the cyclone's 72-h evolu- 
tion is signficantly better than the NMC NOSAT prediction. At 12 h (Fig. 1), 
the error of the cyclone's position is reduced by 51 percent and the cyclone's 
central pressure error is reduced by 66 percent. At 24 h (Fig. 2), no cyclone 
is evident in the NMC prognosis. At 36 h (Fig. 3), the differences between 
the two sea level pressure forecasts are smaller. However, a more intense 
cyclonic circulation is predicted by the GLAS SAT, and the greater amplitude 
of the associated thickness pattern in the SAT forecast is also in better 
agreement with the analysis. At 48 h (Fig. 4), just prior to the outbreak 
of severe thunderstorms and tornadoes, major differences between the GLAS 
SAT and NMC NOSAT prognosis are evident. The structure, intensity, and 
position of the SAT cyclone is substantially better than in the NMC NOSAT. 
The position error has been reduced by 84 percent; the reduction of central 
pressure error is 45 percent. This trend continues to 60 h (Fig. 5), where 
the position error is reduced by 88 percent and the central pressure error 
by 46 percent. At 72 h (Fig. 6), there is a 34 percent reduction in position 
error and 75 percent reduction in central pressure error. 

To investigate the effect of these differences further, the Com- 
puterized Severe Storm Model (CSSM), developed by Atlas (1978), has been 
applied to both the GLAS SAT and NMC NOSAT predictions. The CSSM consists 
of an objective procedure for computing and combining specific measures of 
instability and destabilization to yield a prediction of high-moderate, low, 
or neglible potential for severe local storm development. The SAT prediction 
of high potential accounted for 70 percent of the severe thunderstorm and 
tornado occurrences while the NMC NOSAT CSSM did not predict any areas of 



high potential. The differences between the two CSSM forecasts were due 
primarily to a substantially improved prediction of differential equivalent 
potential temperature advection and convective instability in the SAT system. 

Finally, differences between the GLAS SAT and NMC NOSAT initial 
states have been examined by means of horizontal maps and vertical cross 
sections of a variety of primary variables and derived quantities. Fig. 7 
presents cross sections of absolute vorticity in relation to the initial 
position of the surface low. In the NMC NOSAT initial state, the upper 
level vorticity maximum is located directly above the surface low, whereas 
in the GLAS SAT there is a substantial slope of the vorticity maximum with 
height. The surface low in the SAT analysis is baroclinical ly unstable 
and intensifies as it moves inland during the first 12 h of the forecast, 
whereas the same low fills in the NMC prediction for the same period. 

Forecasts from 0000 GMT 11 February 1976 

Fig. 8 depicts the NMC 7-level NOSAT and GLAS SAT 72-h sea level 
pressure (solid lines) and 1000-500 mb thickness (dashed lines) forecasts 
and corresponding verification. The 24- and 48-h predictions (not shown) 
are contained in Atlas (1979). 

At 24 h, there has been a 66 percent reduction of position error 
and an 87 percent reduction of central pressure error of the cyclone located 
in southwest Canada. Similarly, at 48 h, there has been a 76 percent reduc- 
tion of position error and a 25 percent reduction of central pressure error 
of this cyclone, now centered near the Great Lakes. The position and orienta- 
tion of the pressure trough extending southwestward from this cyclone, and 
the position of the anticyclone located off the east coast of the U.S. are 
better forecast by the GLAS SAT. In addition, a spurious cyclone, forecast 
by the NMC NOSAT to be over southwest Canada and the northwestern U.S., does 
not appear in the GLAS SAT. At 72 h, the GLAS SAT is significantly better 
than the NMC NOSAT in its prediction of the cyclonic circulation over south- 
eastern Canada and the northeastern U.S., the anticyclone over the Midwest, 
the pressure trough extending southeastward from southwest Canada, and the 
shallow pressure trough! ng off the west coast of the U.S. The NMC NOSAT is 
slightly better in its prediction of the cyclone located off the northwest 
coast of the U.S. and also has a 26 percent smaller position error of the 
cyclone now located in eastern Canada. However, the central pressure v error 
of this cyclone has been reduced by 53 percent in the GLAS SAT. ^ 

Two cases where NMC's operational model in 1976 had serious diffi- 
culties in foreasting for the U.S. have been examined. The GLAS model, when 
utilizing initial conditions which included satellite sounding data, was able 
to significantly improve upon the current operational model's predictions. 
Satellite temperature soundings appear to have potential to correct gross 
analysis errors which infrequently occur in the operational NMC system in 
data sparse regions. 
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Fig. 7b. As in Fig. 7a except for GLAS SAT. 
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SUPPLEMENTARY NOTES ON THE DEVELOPMENT AND VERIFICATION 
OF THE AUTOMATED FORECASTING METHOD (AFM) 

(R. Atlas and R. Rosenberg) 


Prior to beginning assimilation experiments with data from the 
January-March 1976 Data Systems Test (DS'’6), research was initiated to de- 
velop improved measures for assessing tne impact of remote sounding data 
on weather forecasting. Previously, numerical measures of forecast accuracy 
such as the RMS error or SI score had been the main tool used for this 
purpose. However, impacts in RMS or SI are often difficult to interpret 
in terms of impact on local weather forecasting (Druyan et al. , 1975; 
Houghton and Irvine, 1976; Sanders, 1973). This is due to the large-scale 
averaging inherent in this type of verification procedure as well as the 
methodology by which weather forecasts are actually produced. Local weather 
forecasting is not based solely on the one or two fields for which numerical 
scores are routinely computed, rather a variety of prognostic primary varia- 
bles and derived quantities are considered in combination (George, 1960; 
Naval Education and Training Command, 1974). Quasi -geostrophic theory forms 
the underlying framework for much of subjective weather forecasting through 
the omega and geopotential tendency equations (Holtan, 1979). 

Recognizing the need for evaluating the impact of satellite data 
on local weather forecasting, a subjective forecast experiment was conducted 
and, utilizing the results of this experiment, the AFM was developed. In the 
subjective forecast experiment, a series of numerical forecasts from initial 
conditions, which included satellite temperature soundings (referred to as 
SAT) and initial conditions which excluded satellite soundings (referred 
to as NOSAT), were generated. For each forecast case, a set of SAT and 
NOSAT prognostic packages (consisting of most of the charts available at 
local weather stations) was delivered to two experienced forecasters. 
These forecasters would then compare prognostic charts, note changes in 
meteorological variables, and then issue a weather forecast at a few pre- 
selected cities scattered throughout the United States. Each forecast was 
based completely on the subjective interpretation of the SAT or NOSAT 
prognostic charts. This experiment revealed that the impact of satellite 
sounding data was large enough to influence local weather forecasts (Atlas 
and Sakai, 1978) but was limited in utility because of the large amount 
of time required to produce weather forecasts at even a small number of 
cities and because of the subjectivity involved. 


The attempt was then made to simulate the forecaster's interpre- 
tation of the complete prognostic package. Objective algorithms were 
developed which correlated well with the subjective judgment of the magni- 
tude and importance of each of the prognostic quantities that were utilized. 
The objective quantities were then combined sequentially in approximately 
the same order as had been employed by the forecasters to form the AFM. 
Minor modifications to the AFM were then incorporated to minimize bias and 
to account for local topographical differences. The AFM is thus based 
primarily on subjective forecasting experience, which is to a large extent 
in agreement with quasi -geostrophic theory. 
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As mentioned in our paper, the AFM was verified in terms of its 
accuracy when applied to a perfect prognoses (i.e., an analysis) and its 
agreement with experienced forecasters, prior to being applied to model 
prognoses. It was found to be 93 percent accurate in determining whether 
or not precipitation would occur, and agreed with the subjective judgments 
of three forecasters (two of which did not take place in the original experi- 
ment on which the AFM was based) 94 percent of the time on the average. The 
AFM was then utilized to evaluate the impact of satellite sounding data 
assimilation on local precipitation forecasts at 128 cities in the United 
States as reported by Ghi 1 et al . (1979). 

Recently, GLAS model forecasts of measurable precipitation and AFM 
forecasts of precipitation occurrence have been compared. Ten different 
72-h forecasts from February 1976 have been verified at 12-h intervals at 
128 cities, for a total of 7,680 precipitation forecast comparisons. Of 
these, the AFM and model forecasts differed 29 percent of the time. Where 
they were different, the AFM was more accurate in its determination of pre- 
cipitation occurrence 67 percent of the time. The greater accuracy of the 
AFM is related to the deficiencies of the model vertical motion predictions, 
as described by Ghi 1 et al . (1979) and to the inclusion of routines to 
account for frontal precipitation in the AFM. 

To determine the effectiveness of the AFM's combination of param- 
eters, we have compared its accuracy when applied to perfect prognoses with 
simple forecast schemes in which the prediction is based on a single parameter. 
Table 1 presents the results of this evaluation for the AFM versus selected 
cutoff values of 500 mb vorticity advection, 850 mb temperature advection, 
or 850 mb dewpoint depression. Comparisons with observed precipitation 
events were made every 12 h for the entire month of February 1976 for the 
same cities that were utilized in our original verification of the AFM. 
Statistics verified included percentual correctness, bias, prefigurance, 
post agreement, and threat score; for a definition of these statistics see, 
for instance, Ghi 1 et al I . (1979) and Tracton and Stackpole (in Murphy and 
Williamson, 1976, p. 736). 

The table shows the relationship of each of the single predictors 
to precipitation occurrence, and the relative utility of each prediction 
scheme. From the prefigurances, we see that 70 percent of all precipitation 
events occurred in conjunction with positive vorticity advection, as compared 
to 49 percent occurring with warm advection and 57 percent with dewpoint 
depressions less than 5°C. However, the forecasts based on weak vorticity 
or thermal advection have higher bias and corresponding lower values of 
percentual correctness, post agreement, and threat score than those based on 
dewpoint depressions less than 5°C. The AFM, with near-perfect bias, accounts 
for 76 percent of all precipitation events and verifies 76 percent of all 
positive precipitation forecasts. This results in a 10 percent increase in 
the total number of correct forecasts, and an increase in threat score of 
.23, relative to the best single parameter prediction. We conclude that: 

a) vorticity advection is a useful predictor of precipitation occurrence, 

b) dewpoint depression is the most accurate single predictor, and c) the 
combination of parameters in the AFM yields a substantial increase in 
accuracy over any single predictor. 
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Table 1. Verification of AFM and single parameter predictions 
of precipitation occurrence for February 1976. 


Forecast 

scheme 

Percent 

correct 

i <as 

Pre- 

figurance 

Post- 

agreement 

Threat 

score 

Precipitation occurs 
when 500 mb vorticity 
advection: 

> 0 x 10" 9 sec" 2 

. 

55 

3.03 

.70 

.23 

.21 

> .5 x 10" 9 sec" 2 

79 

.92 

.35 

.37 

.22 

> 1 x 10- 9 sec -2 

83 

.30 

.17 

.55 

.15 

Precipitation occurs 
when 850 mb temperature 
advection: 

> 0 x 10" 9 sec" 2 

49 

2.91 

.49 

.17 

.14 

> .5 x 10“ 9 sec _2 

74 

.99 

.24 

.24 

.14 

> 1 x 10" 9 sec“ 2 

79 

.39 

.10 

.25 

.08 

Precipitation occurs 
when 850 mb dewpoint 
depression: 

< 5°C 

83 

1.12 

.57 

.51 

,37 

< 2.5°C 

83 

.16 

.16 

,63 

.10 

. . . 

AFM 

93 

.99 

.76 

.76 

.60 
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In summary, the AFM was developed in response to the need for 
improved automated measures for assessing forecasting impact. It simulates 
subjective interpretation of an ensemble of prognostic and diagnostic quanti- 
ties. The AFM has been verified from the standpoint of its accuracy when 
applied to analyses and of its agreement with experienced forecasters. Its 
actual prognostic accuracy was tested for a relatively large number of cases; 
these cases, however, were drawn from a limited number of model forecasts. 
The AFM appears to be comparable in performance to the average skilled weather 
forecaster's interpretations of large-scale numerical guidance and is thus 
a useful tool for evaluating the accuracy and utility of large-scale numerical 
forecasts. 
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FORECAST SKILL AS A FUNCTION OF THE DATA ASSIMILATED 
AND METHOD OF OBJECTIVE ANALYSIS 

(W. Baker, D. Han, and G. Chatters) 


The current skill of numerical weather forecasts degrades rapidly 
in the first 48 to 72 h due to inaccurate initial conditions as well as model 
deficiencies. Recent estimates of the sources of forecast error (e.g. , 
Miyakoda, 1975; Robert, 1976; Somerville, 1976) indicate that inaccurate 
initial data may contribute substantially to the rapid error growth. Errors 
in the specification of the initial conditions are due to both inadequacies 
in the observing systems and data analysis techniques. With the aim of 
extending the limits of predictability, we are evaluating the usefulness of 
new observing systems and the accuracy of different methods of data analysis. 

We have tested two different approaches to the objective analysis 
of the observational data as outlined below: 

Objective Analysis Scheme A 

§ Temperature analysis at the mandatory pressure levels 

• Cressman (1959) weighting with fixed scanning radii 

t No horizontal consistency checks 

• Hydrostatic consistency is not imposed 

• A 5-point smoother is applied to the surface pressure after each 
analysis 

• 12-h assimilation cycle with a +3-hour window used on all 
conventional data. Time-continuous assimilation of satellite 
sounding data (Ghi 1 et al . , 1979). 

Objective Analysis Scheme B 

• Geopotential height analysis at the mandatory pressure levels 

• Cressman (1959) weighting modified such that the scanning radius 
is a function of the data density (Stephens and Stitt, 1970) 

• Horizontal consistency checked after each scan for all data 

• Vertical consistency maintained 

• A smoother (Shuman, 1957) is applied to all fields after each 
analysis 

• 6-h assimilation cycle with a +3-h window used on all data. 
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As may be seen, Scheme B analyzes the thickness of the vertical column, 
maintains horizontal and vertical consistency, determines the radius of 
influence about each model grid point as a function of the data density, 
and utilizes a 6-h assimilation cycle. Scheme A, on the other hand, 
analyzes temperature at each level, does not check for data consistency, 
has a fixed radius of influence for all model grid points, and assimilates 
the conventional data every 12 h with a time-continuous assimilation of 
satellite sounding data (Ghil et al . , 1979). 

Assimilation/Forecast Experiments 

To examine the differences in the objective analysis methods and 
the usefulness *of various data for numerical weather prediction, a series 
of assimilation/forecast experiments has been run, using analysis Schemes 
A and B with different combinations of observational data as summarized in 
Table 1. The first four experiments (1405, 1444, 8449, 8939) utilize the 
2.5° x 3° version of the 2nd-order GLAS GCM (Somerville et al I . , 1974) as 
the assimilation and forecast model. In Experiment 1425, the data are 
assimilated with the 2.5° x 3° GLAS model, and the forecast is generated 
on the 4° x 5° 4th-order GLAS model (Kal nay-Rivas et al . , 1977). The 
2.5° x 3° GLAS model is utilized as the forecast model in Experiment 1298 
for the initial conditions provided by the National Meteorological Center 
(NMC) at 0000 GMT 19 February 1976 from their operational cycle with the 
9-level assimilation model (Stackpole, 1976) using the Flattery (1971) 
analysis. 


Table 1. Description of assimilation/forecast experiments 
for the 0000 GMT 19 February 1976 case. 


' 

1 

Description 

Assimilation 

model 

Objective 

analysis 

scheme 

Forecast 

model 

Satellite 

data 

Aircraft 

temp. 

data 

1405 

NOSAT 

2.5° x 3° 

B 

2.5° x 3° 

No 

No 

1444 

NOSAT 

2.5° x 3° 

B 

2.5° x 3° 

No 

Yes 

8449 

NOSAT 

2.5° x 3° 

A 

2.5° x 3° 

No 

Yes 

8939 

SAT 

2.5° x 3° 

A 

2.5° x 3° 

Yes 

Yes 

1425 

SAT 

2.5° x 3° 

A 

4° x 5° 

Yes 

Yes 

1298 

VTPR+BOGUS 

NMC 

Hough 

2.5° x 3° 

Yes 

No 
















In Experiment 8449, Scheme A objectively analyzed rawinsonde, surface, and air- 
craft wind and temperature data every IE h from 0000 GMT 29 January 1976 
through 0000 GMT 19 February 1976. For Experiment 1444, Scheme B analyzed 
pilot balloon wind data 1 in addition to those data analyzed in Experiment 8449. 

Experiment 1444 began from the initial conditions provided by Experi- 
ment 8449 at 0000 GMT 17 February 1976 with the observational data analyzed 
every 6 h through 0000 GMT 19 February 1976. In Experiment 1405, Scheme B 

analyzed the same observational data as in Experiment 1444 except that aircraft 
temperature data were excluded in Experiment 1405 with the initial conditions 
provided by the NMC operational analysis at 0000 GMT 17 February 1976. For 

Experiment 8939, Scheme A analyzed the same data as in Experiment 8449. In 

addition, satellite sounding data obtained from the NOAA-4 Vertical Temperature 
Profile Radiometer (VTPR) and from the High-Resolution Infrared Spectrometer 
(HIRS) and Scanning Microwave Spectrometer (SCAMS) on NIMBUS 6 were also 
assimilated. 

Evaluation of the Analysis Methods and Observing Systems 

To illustrate the differences in the analysis schemes and the 

importance of various data assimilated, the 500 mb geopotential height and 
wind direction are shown in Figs. 1-3. The differences in the height analysis 
from the final analysis of the assimilation among the various experiments 
should be noted in the enclosed region off the west coast of North America. 
In particular, the amplitude of the short-wave trough is greater in Experi- 
ments 1444 (Fig. 1), 8449 (Fig. 2), and 8939 (Fig. 2) with aircraft temperature 
data than in Experiment 1405 (Fig. 1) where it was excluded. The NMC analysis 
in Experiment 1298 (Fig. 3) also excluded aircraft temperature data, but com- 
parison of Experiment 1298 with the other experiments is difficult because 
the data were assimilated with a different model. Bogus data were also used 
in the NMC analysis in addition to other differences. 

The experiments which utilized Scheme B (1405 and 1444) are less 
noisy than these in which Scheme A was used (Experiments 8449 and 8939). This 
is due in part to the application of a smoother (Shuman, 1957) on the height 
and wind fields in Scheme B, whereas those fields were not smoothed in Scheme A. 
The lack of vertical consistency in Scheme A also contributes to the noise. 

The SI scores for a series of 72-h forecasts generated from the 
initial conditions presented in Fig. 1 are summarized in Table 2. The follow- 
ing points are noted: 

1. Of the five forecasts from initial conditions generated on the 
2.5° x 3° GLAS model (Experiments 1405, 1444, 8449, 8939, and 1425), the 
two which included satellite sounding data and aircraft temperature data 
(Experiments 8939 and 1425) are the most skillful over North America. 


1 For the purposes of our discussions here, the effects of the pilot balloon 
data may be considered negligible. 


Scheme B NOS AT w/o Aircraft Temperatures Exp. 1405 
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A number of points may be made in regard to the diagrams in Figs. 4 
and 5. The long wave error fields shown on the left are quite similar in 
Experiments 1444 (Fig. 4), 8449 (Fig. 5), and 8939 (Fig. 5) during the first 
48 h over Europe (stippled region from 10W to 40E longitude). Aircraft 
temperature data are common to all three experiments. Experiment 1405 (Fig. 4), 
which had neither aircraft temperatures nor satellite sounding data, has sub- 
stantially less error over the same region during the first 48 h of the fore- 
cast. Comparing Experiments 8449 and 8939 over North America (stippled region 
from 130W to 75W longitude), there is no apparent improvement in the forecast 
of the long waves with the addition of satellite sounding data (Experiment 8939). 
However, when Experiments 8449 and 8939 are compared over the same region for 
the baroclinic waves, shown on the right in Fig. 4 and 5, a significant 
improvement is noted, particularly in the last 12 h. This pattern also holds 
for the short waves (k - 8-12) not shown here. 

It is also worth noting that both Experiments 1405 and 8939 are 
generally free of the rapid error growth in the baroclinic waves, character- 
izing the other experiments in the last 12 h of the forecast. It seems likely 
that, when aircraft temperatures and satellite sounding data are assimilated 
together (e.g., Experiment 8939), systematic errors are controlled in the 
barocli nic-wave range and smaller. The rapid error growth seen in Experi- 
ments 1444 and 8449 in waves 5-7 seems consistent with the findings of 
Baumhefner and Julian (1972). In a later study, Baumhefner and Julian (1975) 
found a rapid error growth from initial conditions obtained using simulated 
cloudy retrievals with systematic errors. 
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2. The same pair of experiments, however, are less skillful over Europe 
than Experiment 1405, in which satellite sounding data and aircraft temperature 
data are not assimilated. 

3. Differences in the forecasts with the 4° x 5° 4th-order model and the 
2.5° x 3° 2nd-order model from the same initial conditions (Experiments 1425 
and 8939) are smaller than the differences in the forecasts from different 
initial conditions with the same model (e.g., Experiments 1405 and 8939). 


Table 2. Sensitivity of forecast skill to the data assimilated and method 
of objective analysis for the 0000 GMT 19 February 1976 case. 



Sea level pressure (SI ) 

500 mb i (SI ) 


48h 

N.Am. 

72h 

N.Am. 

m 

■ 

48h 

N.Am. 

72h 

N.Am. 

H 

m 

1405 

76.8 

73.1 

48.8 

67.9 

41.0 

38.3 

<♦6.1 

50.5 

1444 

70.6 

72.2 

53.1 

70.9 

. 

42.5 

45.0 

52.7 

51 .3 

8449 

69.5 

77.5 

52.0 

80.1 

39.5 

44.8 

53.2 

52.2 

8939 

56.4 

54.3 

54.3 

81.6 

33.2 

37.1 

54.5 

62.3 

1425 

61.9 

56.7 

57.9 

82.8 

31.6 

35.1 

48.9 

63.2 

1298 

58.6 

58.8 

48.4 

68.7 

31.2 

35.3 

46.8 

47.2 


Analysis of Forecast Error by Wavenumber Decomposition 

To further examine the effects of the objective analysis schemes 
and data assimilated on the forecast skill, Figs. 4 and 5 show the error 
in the long waves (wavenumbers 1-4) and baroclinic waves (wavenumbers 5-7) 
as a function of time at 40N for 500 mb. The forecast error fields presented 
are relative to the NMC operational analysis and are, therefore, not a measure 
of the "true" observational error. Nevertheless, they provide useful insight 
into the effects of the data and the analysis schemes on the forecast. 
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Longitude-time plot of forecast errors at 40N at 500 mb for 
the 0000 GMT 19 February 1976 case for Experiments 1405, 
1444, and 8449. Errors for waves 1-4 are shown on the left 
and those for waves 5-7 on the right. The contour interval 
is 20 m. 
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Fig. 5. As in Fig. 4 except for Experiments 8939 and 1298 




ASSIMILATION AND FORECAST EXPERIMENTS WITH THE GLAS 
4TH-0RDER GLOBAL ATMOSPHERIC MODEL 

(E. Kalnay-Rivas, W. Baker, and D. Han) 


The new 4th-order GLAS global atmospheric model (Kalnay-Rivas and 
Hoitsma, 1979) is based on an energy-conserving scheme, with all horizontal 
differences computed with 4th-order accuracy. A 16th-order Shapiro filter is 
applied every 2 h to the fields of sea level pressure and potential tempera- 
ture, which are rather insensitive to the presence of topography. In this 
scheme, waves with wavelengths longer than four grid lengths are resolved 
accurately without damping. Waves shorter than four grid lengths, which 
would otherwise be grossly misrepresented by the finite differences, are 
filtered out while they are still infinitesimal. These shortest waves, 
which in the new 4th-order model are considered to be part of the subgrid 
scale, act as a "sponge layer" in the spectral domain, allowing energy to 
trickle down from longer waves and avoiding the accumulation of energy which 
would otherwise occur at the short wave cutoff. Shallow water experiments 
indicate that such a scheme in practice conserves potential enstrophy to a 
high degree of accuracy. 

The new model has similar computation time and memory requirements 
as the present 2nd-order GLAS model with the same 4° latitude, 5° longitude, 
and 9 vertical -level resolution (Somerville et al . , 1974). However, the 4th- 
order forecasting skill is nearly as skillful as the 2.5° lat. by 3° long. 
2nd-order GLAS model. Figs. 1 through 4 show a comparison of the 3-day fore- 
cast from 0000 GMT 19 February 1976 during the DST-6 period made with the low 
resolution 2nd- and 4th-order models, and with the high resolution 2nd-order 
model. A less skillful forecast from 0000 GMT 1 February 1976 is shown in 
Figs. 5 through 8. The relative differences between the models are typical 
of other integrations. The 4th-order forecast for the 19 February 1976 case 
(Fig. 3) was continued for 2 days (Fig. 9), still showing considerable skill. 
Another advantage of this model is that vertical motion is clearly associated 
with only the large-scale field, even over regions with topography (Fig. 10). 
In the 2nd-order models, small-scale noise dominates the instantaneous field 
of vertical motion, especially over topography. 

The analysis and assimilation "Scheme B" (Baker et al . , 1980) has 
been adapted to the 4th-order model . The scheme has been tested with Level Il-a 
(operational) data corresponding to the period 6-9 January 1979 during the 
FGGE SOP-1 . Fig. 11 shows the 300 mb initial analysis at 9 January 1979 
obtained using the 4th-order assimilation model and used as initial conditions 
for the 4-day prediction of Fig. 13. In the eastern hemisphere, the flow is 
rather zonal. After 4 days (13 January), a deep trough developed over Central 
Europe and an intense blocking high, with a low to the south located at about 
60E (Fig. 12). The basic characteristics of this development are well pre- 
dicted by the 4th-order model (Fig. 13). A similar forecast, made from the 
NMC analysis with the high resolution 2nd-order GLAS model, predicts a less 
intense development (Fig. 14). 
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Assimilation and forecast experiments with Level I I-b data using 
both the 2nd~order high resolution and the 4th-order low resolution GLAS models 
are presently being conducted. 
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Fig. 4. GLAS 2nd-order, 2.5° x 3°, 3-day forecast at 0000 GMT 22 February 1976. 






























GEOPOTENTIAL HEIGHT AT 300 



for 0000 GMT 9 January 1979 
















































GEOPOTENTIAL HEIGHT AT 300 MB 



i ; 


rr ^ 

■***<*' s 
caeo- *- 


’r ■ rt ?■ 
t J f 
T $\T 
t ‘ 7" 1 


Sw 

p\'r r 

■t \t \ 

\ K a 

y joV 

k rA" 

2 y\j 

2 \ A 


U^; 3 v» /»' / /' 'pL+V 

Jv F if 

\ A- 

t il-v 

•v <* ySrrfri 

At’v 7 

i \ V 

* \\V 

•yy '•/. y&r' 

7 J r J r , 

/ r r A- 


' s' 


£ 228 ; 


/ 7 r -cot. 
/ W -* 


■t r r. ^ -s 

‘t f V* *•%. 


• * < 
** *✓ 
. *y "/ 

< ✓ 

^ < *" 
* •Vj ^ 

: >< ^ 


- i 


GLAS 2nd-order, 2.5° x 3°, 4-day forecast at 0000 GMT 9 January 1979 











AN IMPLICIT SCHEME FOR METEOROLOGY 
(E. Isaacson and D. Marchesin) 


A compact implicit scheme with accuracy of 4th order in space and 
2nd order in time is developed. The size of the time step is determined by 
the speed of the significant waves in the flow, rather than by the sound 
speed. For a shallow water atmospheric model over a spherical earth, the 
implicit scheme is shown to be more efficient than an explicit scheme, for 
fine grids. 

Meteorological data is becoming available more densely over the earth 
through satellite and other soundings. Numerical methods that can use this 
information to make accurate and longer range forecasts are being intensively 
developed (e.g., Kalnay-Rivas and Hoitsma, 1979). Implicit schemes have the 
advantage that they do not require that the time step be decreased proportion- 
ately when the spatial interval is reduced. That is, the time step of an 
implicit method is chosen to be small enough to resolve time change of the 
motion and not the "sound" speed of the system. An efficient implicit scheme 
for dealing with the horizontal (longitude, latitude) shallow water atmospheric 
model on a spherical earth, with orography, has been tested. For each time 
step, the implicit scheme required slightly more than four times the number of 
arithmetic operations that were needed for each time step of a 4th-order explicit 
scheme developed by Williamson and B/owning (1973). With a spatial resolution 
of (360/128)°, the implicit scheme required a 15-minute time step for a rapidly 
changing flow pattern, whereas the explicit scheme required a 3-minute time 
step, because of the "sound" speed. That is, for this grid of slightly less 
than 3°, the implicit scheme used 0.8 as much computer time as did the explicit 
scheme for a forecast. If the grid size were to be halved, the implicit scheme 
would use 0.4 as much computer time. Thus, for practical application to the 
primitive forecast equations, it is expected that the implicit method will be 
more efficient than the explicit method. 

The shallow water equations can be written in divergence form by 
introducing the momentum variables U = hu and V = hv, where h is the depth 
of the atmospheric layer and u,v the relative tangential velocity components 
from west to east and from south to north, respectively. That is, the equa- 
tions may be written as 


W t + F(W) X + G(W) 0 + E(W) = 0 (1) 


where 



F, G, and E are vector functions of the components of W, and X and 6 are 
longitude and latitude, respectively (see Stoker and Isaacson, 1975, p. 40, 
with R = 0, for the momentum equations in divergence form, without orography). 
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By integrating Eq. (1) from nAt to (n+l)At by the trapezoidal rule, 
one obtains implicit equations involving F, G, and E at Wh+l and at W rl . The 
vector functions may be linearized about W n , which leads to the implicit 
system 


[I + At (JL A” + jL B" + D n )]AW n = R" + 0(At 3 ) 
2 9X 36 


( 2 ) 


where 

AW n » W n+ > - W n , 


A, B, D are the Jacobian matrices of F, G, E, respectively, relative 
to the components of W; R n represents quantities known at time nAt; and the 
partial derivatives are to be applied through AW 11 . 

In any difference approximation of Eq. (2), a two-dimensional system 
of linear equations is found for AW n . But, it would be a formidable task to 
solve this system directly. Instead, by following the formulation of Beam and 
Warming (1976, 1977), the difference operator is factored into the product of 
one-dimensional operators. 

Cl + At (3_ + K n )] [I + ^ (|_ B" + Ln)]Awn - R n + 0(At 3 ). (3) 


where 


K n + L n = D n . 


Experiments were performed with various methods of factoring Eq. (2). 
The following observations were made: 

(a) In some of these methods, the block tri-diagonal, periodic ma- 
trices that result by replacing the derivatives by a 4th-order rational Pade 
approximation have a structure that permits the 3 by 3 blocks to be treated 
as 2 by 2 blocks along with a scalar system. 

(b) Furthermore, although the block systems are periodic, the detec- 
tion of underflows in the bottom row and rightmost column of their LU decomposi- 
tion makes it possible to perform the Gaussian elimination process in only 
about 10 percent more time than would be required for the Dirichlet data block 
system. 

These two observations are primarily responsible for the fact that 
the implicit scheme required only slightly more than four times the number 
of arithmetic operations as are required by the explicit scheme. A 4th-order 
Shapiro filter is applied after each of the spatially differenced factors is 
"inverted." , Si nee this filter is applied to the "change," rather than to the 
function W n+I , no apparent dissipation of the solution is detected (see 
Kalnay-Rivas and Hoitsma, 1979, for discussion of such dissipative effects). 
A recent paper by Navon and Riphagen (1979) carries out a similar implicit 
procedure for the shallow water equations over a rotating rectangle, with 
periodic boundary conditions. 
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Test Calculations 

The implicit scheme was tested with three kinds of initial data: the 

steady zonal flow, the nearly geostrophically balanced initial data derived 
from prescribing three subtropical highs symmetrically placed in northern and 
southern hemispheres, and a four-mode Haurwitz wave (see Phillips, 1959; 
Sadourny, 1975). 

The flow resulting from the initial hi$hs quickly developed lows 
(one for each high and bordering it on the east side) and the pattern moved to 
the west. This rapidly changing flow was a more severe test for the implicit 
scheme. The accuracy of the solution was estimated by comparison with the 
solution produced by an explicit scheme. It was found that the solution of 
the implicit scheme at a time step of 15 minutes agreed quite well (by overlay- 
ing contour maps) with the solution of the explicit scheme. Fig. 1 shows the 
contours of the initial height field with arrows indicating the initial velocity 
field, in a stereographic projection of the northern hemisphere. Fig. 2 shows 
the contours of the height field and the velocity arrows after 1 day as produced 
by the implicit scheme. In these calculations, no mountains were present. 

Summary 


An efficient two-dimensional implicit scheme has been developed for 
the shallow water model over the earth. The scheme conserves total mass, 
energy, and potential enstrophy reasonably well and, when it is required to do 
so, a corrective modification step can be added to ensure proper conservation 
as indicated in Isaacson et al, (1979). Application of the implicit technique 
to the primitive equation model is being developed. Here it is planned to use 
time-splitting to aid in separating the vertical dimension from the tangential 
dimensions. 
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Fig. 1. 


Contour levels of initial height field, from h = 8500 m 
to h = 8585 m. Maximum initial geostrophic velocity is 
about 13 m/s. Arrow size is proportional to speed. 





A LANGEVIN EQUATION FOR LARGE-SCALE ATMOSPHERIC FLOW 
(M. Ghil and R. Balgovind) 


The purpose of this study is to derive on a theoretical basis optimal 
interpolation coefficients for the assimilation of the First GARP Global Experi- 
ment (FGGE) data. In the process, hopefully, a better understanding of the 
stochastic-dynamic behavior of large-scale motions of the atmosphere will be 
obtained. Previous attempts at modeling the statistics of atmospheric fields 
in a systematic way by multi -dimensional autoregressive processes have been 
made by Thiebaux (1977, and references therein). A general theory of statis- 
tical interpolation on the sphere with application to meteorological problems 
is in the process of being formulated by Wahba (1979, and references therein). 

The short-time dynamics of extratropical large-scale atmospheric 
flow are governed by the conservation equation of potential vorticity: 

(d/dt)(v2-x-2)<j» = 0 . (1) 

In it, d/dt is the material derivative, V 2 the horizontal Laplacian, A the 
Rossby radius of deformation, and ^ the height of the 500 mb pressure surface. 
The radius A, a 2 = <J>/f2, represents a characteristic length scale of quasi- 
geostrophic flow, with <J> a characteristic magnitude of $ and f = 2ft sin<J> the 
Coriolis parameter, ft being the angular velocity of the earth and 0 latitude. 

A stochastic perturbation of this equation on the sphere, for an atmosphere 
near geostrophic equilibrium, is considered: 

(v 2 - *-2)* ■ W , (2) 

where W is uncorrelated white noise representing smaller scale nongeostrophic 
phenomena ; more precisely, E W(v»0) = 0, EW(v+v'.S+ 0 ' )W(v ,0 ) =a 2 S(v 1 ,0 ' ) » , 
with E the expectation operator, v longitude, and a 2 a prescribed variance. 

Numerical solutions of Eq. (2) for different values of the parameter 
are obtained. A comparison is made of Rn (v';0) = En<J)(v+v', 0) <J>(v»0) with 
the corresponding spatial correlations of the atmospheric temperature field, 
obtained from satellite observations during the Data Systems Test, DST-6 
(February 1976), for FGGE (Ghil et al . , 1979). Here En is the finite-sample 
correlation based on N realizations of the equation. 

The major feature of the data (Fig. 2 of Ghil et al . , 1979), namely, 
a general R^(s;0) ^exp(-s/s 0 ( 0) ) behavior, is reproduced by the solutions of 
Eq. (2), with s 0 ( 0 ) ~ f - si n 0. For realistic values of $, the fit between 
solution correlations and data correlations is rather good (Fig. 1). This 
seems to indicate that a single vertical mode strongly dominates the error 
structure of the temperature field. 

The theoretical correlation function for Eq. (2) has also been 
obtained, i.e. , R(s; 0 )=1 imj^Rf^s; 0) . This will be further compared with 
the error structure of other atmospheric fields, namely, geopotential heights 


at 5°0 mb and at 850 mb, and sea level 
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APPLICATIONS OF ESTIMATION THEORY TO NUMERICAL WEATHER PREDICTION 
(M. Ghil, J. Tavantzis, K. Bube, and E. Isaacson) 


Numerical weather prediction (NWP) is an initial-value problem for 
a system of nonlinear partial differential equations (PDEs) in which the 
initial values are known only incompletely and inaccurately. Data at initial 
time are supplemented, however, by observations of the system distributed 
over a time interval preceding it. Estimation theory has been successful in 
approaching such problems for models governed by systems of ordinary differ- 
ential equations (ODEs) and of linear PDEs. An attempt to apply it to NWP is 
being made. 

A model exhibiting many features of large-scale atmospheric flow 
important in NWP is the one governed by the shallow-fluid equations. Cur- 
rently, the estimation problem for a linearized version of these equations is 
being studied (Ghil and Tavantzis, 1978). A discrete Kalman-Bucy (K-B) filter 
is applied to a finite-difference version of the equations, which simulates 
the numerical models used in NWP. The filter has been modified in order to 
el imi nate fast inertia-gravity waves in the solution, as well as reduce the 
error in the estimate of the atmospheric state. The vector of observations 
corresponds to the different atmospheric quantities measured and space-time 
patterns associated with conventional and satellite-borne meteorological 
observing systems. 

The specific character of the equations' dynamics gives rise to the 
necessity of modifying the usual K-B filter. The modification consists in 
eliminating the high-frequency inertia-gravity waves which would otherwise 
be generated by the insertion of observational data. The modified filtering 
procedure developed here combines in an optimal way dynamic initialization 
(i.e. , elimination of fast waves [Baer and Tr i bbi a , 1977; Ghil et al . , 1977; 
Ghil, 1980]) and four-dimensional (space-time) assimilation of observational 
data (Bengtsson, 1975; Ghil et al . , 1979; Bube, 1978), two procedures which 
traditionally have been carried out separately in NWP. 

Optimal error reduction in the initial-state estimate is studied 
for the different observation patterns. Comparisons between the modified 
filter and the standard K-B filter have been made. Also the filter discussed 
will be compared with operational methods of four-dimensional data assimilation 
currently used by meteorological services. 
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INTEGRATED AND SPECTRAL ENERGETICS OF THE GLAS GCM 
(J. Tenenbaum) 


Analysis of the generation, transport, and dissipation of energy 
represents one of the primary methods of studying geophysical fluid flows. 
Lorenz proposed subdividing atmospheric energy into kinetic and available 
potential forms and in turn subdividing these into zonal kinetic (KO) and 
eddy kinetic (KE) categories. Further subdivision in the spectral domain was 
suggested by Saltzman (1970) and followed by similar observational and model 
studies (Tenenbaum, 1976; Baker et al . , 1978). One purpose of this project 
is to examine both the integrated and spectral energetics for two versions of 
the GLAS general circulation model. 

The most direct test of a model over a period of several days is 
the quality of the prediction of standard meteorological variables. One 
knows that this approach should theoretically give an essentially random 
result at a predictability limit on the order of 14 days, and in practice 
does so in about 5 days. 

To extract further information in spite of these limits, one employs 
various space and time averages applied to both meteorological variables and 
energetics quantities. These averaged quantities may show skill for periods 
approaching the predictability limit; indeed, current general circulation 
models are characterized by quite reasonable model "climates," In this pro- 
ject, we will examine to what extent low wave number spectral quantities may 
qualify as "averaged" quantities. 

Central results of this project concern the interplay of the eddy 
kinetic energy and its spectral subcomponents. Figs. 1 and 2 show the 
observational and model values of KO and KE for the first 2 weeks of each 
time period for two sets of model runs. Figs, la and lb refer to a GLAS 
model with no filtering of diagnostic variables and a radiation routine 
comparable to the GISS version of the model (Somerville et al . , 1974). 
Figs. 2a and 2b refer to a GLAS model with filtering of diagnostic variables 
and a new radiation routine due to Wu. 

For the no-filtering version, the KE agreement for January 1977 
(Fig. la) is quite good. Note that for the KE rise in January 1977 there 
appears to be a 1-day lag. A reasonable interpretation is the clear evidence 
of predictive skill in an averaged quantity out to the predictability limits. 
With the same model for January 1975 (Fig. lb), the observational KE dip is 
not followed downward. 

For the filtered version of the model, this same case (January 1975) 
is well modeled (Fig. 2a). The accurate tracking of the dip in KE is a 
result which is consistent with the filtering out of extraneous small-scale 
energy. A second filtered case, however, shows far less success: for January 

1976 (Fig. 2b), the model seems unable to generate the KE needed to follow 
the atmosphere's rise. 
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Time history of the tropospheric mean kinetic energy, KO, and eddy 
kinetic energy, KE, for the first two weeks of each case using an 
unfiltered model, (a) January 1977 (b) January 1975. Units 10 5 J/m 2 








One uniform pattern does appear in KO for both versions. The 
initial states consist of one low KO state (January 1976), one average KO 
state (January 1975), and one strikingly high KO state (January 1977). The 
latter case is unusual and was sustained at, or near, this high value 
throughout the 2 weeks. Both models, in contrast, for the four cases have 
only one pattern: monotonic growth. Some of the success of January 1977 

is clearly due U the atmosphere happening to match this built-in bias of 
the model, 


Upon studying spatial displays of KO, excess energy is found 
collecting in the model's topmost layer. Since there are no significant 
dissipative mechanisms at that level, the excess energy appears to remain 
trapped there, functioning as a source of the model's total KE, but feeding 
the individual wave numbers erroneously. As shown, a related study in con- 
junction with E. Kalnay-Rivas of GLAS, this effect is an initial condition 
shock which induces major oscillations over periods from 2 h to 4 days. 
Enhanced vertical resolution near the tropopause decreases the amplitude of 
the initial shock. 

In striking contrast to the successful KE predictions for January 
1975 and January 1977 is the behavior of the individual K n , the coefficient 
of kinetic energy fuf wavs number n. One can isolate the one or two dominant 
wave numbers with little difficulty. Fig. 3 shows the time history of Ki and 
Ko for January 1977, In spite of the noisier signal, one sees a much more 
discrepant behavior with gross errors at times when the summed KE remains 
correct. Note particularly the behavior of n - 3. In terms of model develop- 
ment, the fundamental question is "What process or processes in the model are 
producing the difficulties?" 

In terms of Saltzman's spectral energetics analysis, one can imagine 
two possibilities: (1) erroneous transports in or out of the northern hemi- 

sphere troposphere, or (2) erroneous conversions from zonal potential or 
kinetic energy. A third possibility, erroneous transports withi n the northern 
hemisphere troposphere, should not logically affect the results since inte- 
grated quantities are considered. Such erroneous transports may, however, 
be suggestive of other problems. Inevitably, the sorting out of cause and 
effect will probably be difficult. 

For January 19/7, there is evidence of transport difficulties. 
Northward transport starts in the subtropics and feeds first the large K3 
peak at 38N, and then the stratospheric warming at 50N. For January 1975, 
evidence indicates erroneous conversion in the term expressing nonlinear 
redistribution 0* energy to wave number 3 from all other wave numbers 
(Saltzman's L). At the time (2 January 1975) when K3 for observation and 
model are just starting to diverge, one sees a discrepancy centered at 36N 
and 285 mb. No corresponding effect appears before the much stronger n = 3 
divergence of 6 January 1977. No other conversion term shows significant 
causal discrepancies though a fair number of correlations are observed, i.e., 
D3 (the dissipation in wave number 3) showing significant disagreements when 
K3 does. 



Fig. 3. Time history of eddy kinetic energy of wavenumbers 1 and 3 for 
January 1977. Units 10^ J/m^. 
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Some general questions remain concerning the consistency of results 
obtained from Sal tzman's spectral energetics analysis. The equations themselves 
are not in question but rather the quality of results from their application 
to either model or observational data. Saltzman's equations provide a closed 
system with respect to conversions into and out of individual K n » The algebraic 
sum of the conversions integrated over an appropriate space and time period 
should equal the secular change in the corresponding K n . While such a result 
requires assuming the validity of dissipation calculations, D, the lack of 
agreement gives some sense of overall confidence in the results. Because of 
the unavailability of independent data for D, Saltzman had to calculate it 
only as a residual . 

The results are not encouraging. Table 1 presents this type of cal- 
culation for a sampling of model and observational runs. For most cases, the 
sign of the change is incorrect and, in some cases, the individual K n should 
have been driven solidly and unphysically negative. Three possible explana- 
tions seem most likely: errors in Cn, Dn, and in the vertical transports. 

The conversion Cn and, to a lesser extent, some others depend on the vertical 
velocity, W. For the observations, this is a notoriously difficult quantity 
to evaluate. Dn is a parameterized quantity whose physical basis is known in 
broad terms but whose detailed quantitative evaluation requires confidence in 
surface winds and boundary layer effects. Vertical transports are clearly a 
problem for the model and may be a problem for the observations. 

The vertical velocity problem is being examined in detail in a subse- 
quent paper (Tenenbaum and Reddi , 1980, to be published). Preliminary results 
show that alternate formulations for omega (changing from mass convergence to 
iterative solution of the diagnostic omega equation) can yield major changes 
in long- and synoptic-scale values of Cn and extensive smoothing of the higher 
wave numbers. 

The model's medium-range skill and climatology using various possible 
integrated quantities have been studied. Major conclusions are as follows: 

(1) The model can show medium-range skill in an integrated quantity, 
KE, the northern hemisphere tropospheric eddy kinetic energy, out to the pre- 
dictability limit of 2 weeks. When the atmosphere's K0 matches the model's 
bias towards monotonic growth in K0, the predictions are most successful. 

(2) The model retains a bias in accumulating K0 in the topmost, 
stratospheric level and in the initial monotonic growth of tropospheric K0. 

(3) Low wave number spectral coefficients appear questionable in the 

role of integrated energetics quantities: KE can show agreement while indi- 

vidual Kn are erroneous. An expanded version of Saltzman's spectral analysis 
(including transports) does provide a useful tool in showing cause-and-effect 
relations for where the model is failing to match significant energy flow in 
the atmosphere. In the near future, more synoptic cases will be studied, 
both randomly chosen, and specially chosen to explore synoptic dependence. 
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Table 1. Conversions to and from K3 for wave number 3 for January 
1975 (5 to 7 days) and January 1977 (3 to 6 days). 


Conversions to and from 

January 1975 
(5 to 7 days) 

January 1977 
(3 to 6 Days) 

K 3 (W/m 2 ) 

Observations 

Model 

Observations 

Model 

+C(P 3 ,K 3 ) = c 3 

-0.3 

0 

-0.5 

0.3 

+C(K n ,K 3 ) - L 3 

0 

-0.3 

0.3 

0 

“03 

0.1 

0.1 

0.25 

0.25 

-C(K 3 ,Ko) = M 

0.2 

0.4 

0.1 

0.3 

Net conversion (W/m 2 ) 

-0.6 

-0.8 

-0.55 

0.25 

X time (sec) 

1.7 x 10 5 

1.7 x 10 5 

2.6 x 10 5 

2.6 x 10 5 

Net conversion (10 5 J/m 2 ) 

-1.0 

-1.4 

-1.3 

0.7 

K 3 (final) 

1.5 

1.5 

2.5 

0.5 

-K 3 (initial) 

0.6 

1.5 

1.5 

1.2 

Net change (10 5 J/m 2 ) 

0.9 

0 

1.0 

-0.7 
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FRACTIONAL CLOUDINESS AND VARIABLE 
CLOUD EMISSIVITY IN THE GLAS GCM 

(G. F. Herman and V. Krishnamurthy) 


The recent analysis of cloud-radiation experiments conducted with the 
GLAS GCM ( cf. Herman et al . , 1980) clearly pointed to the fact that the average 
infrared opacity of the GLAS model atmosphere was too large. This is attributed 
to two factors: First, since all model clouds, including those which are intended 
to represent cirrus clouds in nature, are assumed to have unit emissivity (zero 
transmissivity), the model did not transmit enough long wave radiation to space. 
Secondly, all model clouds, including those intended to represent subgrid scale 
cumulus convection, were assumed in the radiation calculation to occupy the 
entire 4° latitude by 5° longitude grid area. The net result of both of these 
features of the radiation calculation was to maintain an excessively positive 
global radiation balance, and possibly have a detrimental effect on such, 
dynamical processes as the formation of the Siberian anticyclone. 

If the spectra of atmospheric absorbers, namely water vapor and 
cloud droplets, are uncorrelated (cf. Goody, 1964, p. 123), then a particularly 
simple treatment of both variable cloud emissivity and fractional cloudiness 
follows from the so-called multiplicative property . 

Suppose that f is the fraction of a grid square that contains ajn^ 
cloud at an^ level, and also that f is constant with respect to height. Then 
the flux.Tc", at any height, z, in this cloudy region is 

00 

F c (z) = Tc(°»z) ir (B(0) - Bg) - y'TcU.z' ) T , dB( z * ) h? 1 . (1) 

o dz' 


where T c is the frequency-averaged transmittance that includes the presence 
of the cloud droplets, B is the Planck intensity, and B g corresponds to the 
temperature of the ground. 

For cloud-free regions of the atmosphere, the net flux due to gaseous 
absorption, F g , is 


F g (z) = T g (0,z) tt(B( 0) - B g ) - /TgU.Z 1 ) p dB.Cz,;,) dz - , (2) 

o dz ' 

where the mean transmittance, T g , accounts only for gaseous absorption. We 
define the average flux, F, of the partially cloudy region as the weighted 
mean, 


F = fF c + (l-f)F g 
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Then from (1 ) and (2) 


F « it ( B (0) - B g ) [f T c (0,z) + (l-f,'Tg(0,z)] 

/ “ (4) 

* dB(z;W fT c (z,z' ) + (l-f)Tg(z,z*)] . 

0 dz' 

The bracketed term CfT c + (l-f)Tg] represents the average transmittance of 
the partially cloudy regions. 

Using the multiplicative property, the mean transmittance, T, of 
the cloud-plus-gas mixture is 


T(z,z')=T g (z,z’)xT c (z,z') . (5) 

The effective transmittance, T e , of the partially clouded region is 

thus 


T e =T g (l-f+fT c ) . (6) 

Equation (4) becomes simply 

F = tt(B(0) - B g )T e - / dB(z') TeU.z 1 ) . (7) 

o 


Equation (7), as derived, applies only to the situation where cloud 
fraction does not vary with height. However, it is readily generalized to 
account for cloud fractions that may vary from layer to layer. 

__ The transmittance of the cloud may thus be adjusted directly jn the 
term T c . Note that the cloud fraction and cloud transmissivity occur as a 
product; the average opacity of the atmosphere may be adjusted by alternatively 
varying the cloud fraction or the cloud emissivity. 

Table 1 illustrates the effect that variable cloud fraction and 
cloud transmissivity have on fluxes and cooling rates. Here f denotes the 
fractions of a grid square that contains any cloud, and T c the cloud trans- 
missivity. Thus, f=0 corresponds to a completely clear grid, while f =1 . 0 
corresponds to a complete overcast. Similarly T c =0 corresponds to a completely 
transparent cloud, while T c =1.0 denotes a cloud that is perfectly opaque. 
For illustration it is assumed that a single cloud layer is found in layer 5. 
The effect of the cloud on the fluxes is shown in upper section of Table 1. 
As the average opacity increases due to either increasing the fraction or 
decreasing the transmittance, the net flux beneath the cloud layer decreases 
monotonically, which is due to the enhanced downward flux. The flux above 
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the cloud layer also diminishes with increasing opacity because the upward 
radiation is less intense, having originated at comparably colder temperatures. 

The effect on cooling rates is illustrated in the lower section of 
Table 1. Cooling below the cloud is smaller as the opacity increases because 
of the decreased importance of cooling to space. Cooling at and above the 
cloud increases with increasing opacity. These dependencies are illustrated 
in Fig. 1. 


When introduced to the GLAS GCM, the cloud fraction, f, will likely 
be a parameter to which the model's radiation and dynamics are very sensitive. 
Clearly, a major problem will be in determining the most realistic value of 
this parameter. 
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Table 1. Clouds in layer 5. Latitude * 46°, Day * 1 February 


Level 


FI ux(W in - 2) 

Pressure 

(mb) 

Cloudless 
f»0.0,T c *l .0 

f*Q.25,T c *0.0 

or 

f*l .0 ,T C B 0.75 

f=0.5,T c s! 0.0 

or 

f-1 .0,T C =0.5 

f e 0.75,T c *0.Q 

or 

f c *l •°»T C *0.25 

f-1 .0 
T c *0.0 

1 

10 

211 

195 

179 

164 

148 

2 

120 

205 

189 

173 

157 

141 

3 

230 

197 

181 

164 

148 

131 

4 

340 

187 

170 

153 

136 

119 

5 

450 

170 

152 

134 

116 

98 

6 

560 

151 

126 

100 

75 

49 

7 

670 

134 

113 

92 

72 

51 

8 

780 

118 

101 

83 

66 

48 

9 

890 

104 

89 

74 

59 

44 

10 

1000 

90 

77 

64 

51 

39 


Layer 


Cooling Rate (°C/Day) 

Pressure n 
at Center 
of Layer 
(mb) 

Cl ou^less 
f=0.0,Y c =1.0 

f=0.25,T c = 0-0 

or 

f=1.0,T c =0.75 

f=0.5,T c =0.0 

or 

f=l .0,T C =0.5 

f=0.75,T c =0.0 

or 

f=l .0 ,T c =0.25 

f-1 .0 

T c =o.o 


65 

0.4 

0.5 

0.5 

0.5 

0.6 


175 

0.6 

0.6 


0.7 

0.7 


285 

0.8 

0.9 

0.9 

0.9 

0.9 


395 

1.3 

1.4 

1.4 

1.5 

1.6 


505 

1.4 

2.0 

2.6 

3.2 

3.8 


615 

1.3 

1.0 

0.6 

0.2 

-0.2 

7 

725 

1.2 

1.0 

0.7 

0.5 

0.2 

8 

835 

1.1 

0.9 

0.7 

0.5 

0.3 

9 

945 

1.1 

0.9 

. . . 

0.8 

0.6 

0.4 


^f**^* 8 *'** 1 *" ww®**^* 
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ABSORPTION OF SOLAR RADIATION IN THE GLAS GCM: 
THE INFLUENCE OF SOLAR ZENITH ANGLE 

(R. Davies) 


The current parameterization of the bsorption of solar radiation 
in the GLAS GCM is relatively unchanged from the initial scheme developed 
by Lacis and Hansen (1974), referred to hereafter as the LH code. While 
the LH code explicitly includes the effect of the solar zenith angle, 6 0 , 
on the absorption of solar radiation in cloud-free atmospheres, it fails to 
include this effect in the two critical areas of surface absorption and 
cloud scattering. The consequences of such an omission with regard to 
climatic effects and to model development are discussed below. 

Surface Absorption 

Following the observations of many investigators, which have been 
summarized by Paltridge and Platt (1976), it is now recognized that the 
surface albedo of most natural surfaces under clear sky conditions shows a 
pronounced dependence on 0 O at iow solar elevations. If an is the clear 
sky surface albedo for 0 o <40°, this dependence may be modeled as 


a g (0 O ) * *g + 0 - a g )exp{-0.1 (90° - e 0 )} . (1) 


so that a g ->-1 as 6 0 -*90°, but remains close to an for 0 O <40°. While this 
effect may be incorporated into the GCM by a relatively trivial change to 
the code, the consequences to the model -generated planetary albedo and 
surface heating are quite significant. One should note that Q 0 is greater 
than 40° for at least 50 percent of the day, even at equatorial latitudes. 

The most pronounced effect of this change is on the amount of 
absorbed radiation by the earth's surface. Using model -generated values 
for cloud, water vapor, etc., at 1200 GMT 1 January, zonally averaged 
values of surface absorption were produced using the LH code and the new 
code (NC). Fig. 1 shows the difference values (LH results - NC results) 
as a function of latitude. 

The LH code systematically overestimates the surface absorption 
by up to -45 ly/day, with the worst bias at Mgh latitudes. On a global 
basis, the bias is ~30 ly/day for this particular case, representing an 
error of ~5 percent. If this difference were to be applied continuously to 
the snow line, the LH code would melt an extra inch of snow every 2-3 days. 
This error may also be partially responsible for the abnormally high surface 
temperatures produced by the GCM at high latitudes. 
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Cloud Scatteri ng 


The effect of solar zenith angle on the transfer of solar radiation 
in an isolated plane parallel cloud layer with homogeneous intrinsic parameters 
t, a , and g is modeled first. Here x is the optical thickness, a is the single 
scatter albedo and g is the asymmetry factor of the scattering phase function. 
The required output quantities are the amounts of solar energy absorbed, trans- 
mitted and reflected by the cloud layer. Distinction will be made between two 
forms of these quantities; those arising from a direct sola; beam of 'rradiance 
F 0 from zenith angle 0 O , and those arising from unit diffuse irradiance on the 
top surface of the layer. The total problem of multiple cloud layers, cloud- 
free intervening atmosphere and underlying surface may then be solved by a 
suitable combination of the isolated layer results, as discussed below. 

In modeling the transfer of solar radiation in clouds, the LH code 
uses a diffuse version of the two-stream approximation (Sagan and Pollack, 
1967), whereas the new code employs the delta Eddington approximation (Joseph 
et a!. , 1976) for large x and a delta version of the Coakley and Chylek (1975) 
two-stream approximation for small x. Both the delta Eddington and the delta 
two-stream models explicitly include solar zenith angle effects. 

Fig. 2 presents comparative results from the LH and new codes for 
the transmitted energy through a conservative cloud of optical thickness 
x = 30 as a function of Pq = cos 9 0 resulting from an incident direct beam 
with F 0 - 2880 ly/day. Also (,wm are results from the accurate doubling 
model of Hansen (1971). While the delta Eddington model provides excellent 
comparative results to the doubling model, the LH version is seen to system- 
atically underestimate the transmission for large values of u 0 and to over- 
estimate it for low values of y 0 . Due to the nonlinearity of the radiative 
transfer problem, it appears unwise to rely on averaging to cancel the errors, 
as was suggested by Lacis and Hansen (1974). 

Figs. 3 and 4 address another important feature of cloud scattering, 
namely its effect on absorption within an individual cloud layer. The absorbed 
energy is presented as a function of t for a = 0.7 and pq = 1.0 in Fig. 3, and 
as a function of a for x = 3 in Fig. 4. The new code and the doubling model 
give virtually identical results in Fig. 3 and excellent agreement in Fig. 4, 
whereas the LH code continues to show systematic errors in both figures. 

While the LH code failed to adequately account for the directionality 
of the incident solar beam, it gave relatively good values for the transmitted, 
reflected and absorbed flux resulting from diffuse incident irradiance. An 
example of this is shown in Fig. 5, which compares the results for the trans- 
mitted fraction of diffuse irradiance from the LH version and the doubling 
model as a function of x for different a, 

In developing the new code, solutions to the complete multi-layered 
problem were Stained by a combination of the individual layer solutions. 
This was accomplished in two stages. Upwelling and downwelling irradiances 
from each cloud layer and the upwelling irradiance from the surface due to 
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i. 2. Transmitted solar flux as a function of Vo* *^ e 
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Fig, 3. Absorbed solar irradiance as a function of cloud 
optical thickness for F 0 * 2880 ly/day. 




Fig, 5. Relative transmission of solar flux for diffuse 

incidence as a function of cloud optical thickness 



the attenuated solar beam at each level were first calculated. The resulting 
interaction of these irradiances with the complete multilayered system was 
then obtained by applying an advanced adding technique, utilizing the diffuse 
coefficients for reflection and transmission. The solar zenith angle depend- 
ence was thus preserved while still allowing the complete problem to be solved 
in a computationally efficient manner. Preliminary estimates using the new 
code to generate the global results below, for example, indicated an increase 
of less than 10 percent in the running time of the solar radiation portion of 
the GCM. 


On applying the revised treatment to a typical time step of the 
GCM, the main effect on the zonal ly averaged solar radiation was an increase 
in planetary albedo and a decrease in atmospheric absorption. Fig. 6 shows 
an example of the difference between the new and old codes in generating 
total atmospheric absorption. The original treatment of surface albedo was 
used in this comparison, so that the differences resulted entirely from the 
new treatment of radiation in clouds. The values were zonal ly averaged for 
the same case as in Fig. 1. The LH code again systematically overestimates 
absorption for most latitudes, but with a latitude dependence that strongly 
reflects the amount of cloud cover. The maximum difference is thus now in 
equatorial latitudes where there is both significant cloud cover as well as 
substantial insolation. On a global basis, the LH code overestimates the 
atmospheric absorption for this case by ~20 ly/day or ~20 percent,. 

The present study has demonstrated that the absorption of solar 
radiation both at the earth's surface and within the atmosphere is noticeably 
affected by the directionality of the insolation at the top of the atmosphere. 
The combined effect of the above revisions on the zonal ly averaged planetary 
albedo is presented for the 1 January case in Fig. 7, which shows both 
the new and LH results. The increased surface albedo for clear skies and 
the increased cloud albedo for cloudy skies resulted in an increase of the 
planetary albedo at virtually all latitudes. On a global basis, the albedo 
for this particular case rose from the LH value of 0.29 to 0.34. In terms 
of absorbed solar energy, this has an equivalent climatic effect of a decrease 
of 7 percent in the solar constant. 

Further developments along the above lines are planned to more 
accurately account for spectral effects and to more faithfully model situa- 
tions involving fractional cloudiness. 

Acknowledgment 

Financial support under NASA Research Grant NSG 5223 is gratefully 
acknowledged. 

References 

Coakley, J. A., and P. Chylek, 1975: The two-stream approximation in radiative 
transfer: including the angle of the incident radiation. J. Atmos. 
Sci . , 32, 409-418. 

Hansen, J. E., 1971: Multiple scattering of polarized light in planetary 

atmospheres. Part I. The doubling method. J. Atmos. Sci., 28 , 120-125. 


Planetary albedo oiffehence in aimosphcric AnsonrnoN 



-io r 


-so 1 " 
s 

Fig. 6. GCM results for the difference In zonal ly averaged 
atmospheric absorption (LH results - new results). 
Model time is 1200 GMT 1 January. 


-6 O’ 


-30" 


O' 

LATITUDE 


j .... 

3 0' 


60' 


90 

N 



S LATITUDE N 

Fig. 7. Zonal ly averaged planetary albedo produced by 
the LH and new codes in the GCM. Model time is 
1200 GMT 1 January. 

102 




f 


\ 


Refere nces (Continued) 


Joseph, J. H., W. J. Wiscombe, and J. A. Weinman, 1976: The delta-Eddinqton 

approximation for radiative flux transfer. J. Atmos. Sci. . 33, 2452-2459. 

Lacis, A. A. , and J. E. Hansen, 1974: A parameterization for the absorption of 
solar radiation in the earth's atmosphere. J. Atmos. Sci. . 31, 118-133. 


Paltridge, G. W. , and C. M. R. Platt, 1976: Radiative Processes in 

and Climatology , Elsevier Sci. Publ. Co.TWwTorTTT 


Meteorology 


Sagan, C., and J. B. Pollack, 1967: Scattering in the Venus clouds. J. Geoph.ys. 

kgs » f ! c. y Hoy "*♦/ / • ■ 


103 




TESTS OF GROUND HYDROLOGY PARAMETERIZATIQNS 
(0. D. Lin, P. Bock, and 0. J. Alfano) 


A global ground hydrology model (GHM) is being developed for use in 
the GLAS GCM. Parameterization of land surface hydrologic processes at the 
grid scale of the GCM is the central problem in the model development and 
research. Ground temperature and soil moisture are employed in the GCM as 
prognostic variables and thus must be used as basic hydrologic variables in 
the GHM for parameterization of the latent heat flux (evapotranspiration) and 
sensible heat flux at the atmospheric-land Interface. To account for the 
pathways through which the moisture and hwc are transported, a computational 
grid in the GHM is divided horizontally into cuo portions: the bare soil and 
the vegetated surface, and vertically into Wo active layers: the surface 
layer and the root zone, as depicted schematically in Fig. 1. The subgrid 
parameterization is considered to be -a minimum requirement for a physically 
realistic GHM. The purpose of this report is to present the test results of 
the GHM by using the field-observed data that were obtained from two experi- 
mental sites on bare soils and one on winter wheat crop. The data consist of: 

(1) soil moisture, ground temperature, and evaporation of bare 
soil from the U. S. Water Conservation Laboratory, Phoenix, Arizona (Jackson, 
1973; Jackson e t al . , 1975). 

(2) soil moisture and ground temperature of bare soil from the 
Southern Piedmont Conservation Research Center, Watkinsville, Georgia (Bruce 
et al . , 1977). 


(3) soil moisture, canopy temperature, and evapotranspiration of 
winter wheat crop from U. S. Water Conservation Laboratory, Phoenix, Arizona 
(Jackson et al . , 1977). 

The formulations that are involved in the one-dimensional transport 
of moisture and heat in the tests of a soi 1 -atmosphere or soil -plant-atmosphere 
column are summarized as follows. A detailed presentation of the GHM is given 
in Lin et al . (1979). The ground temperature is calculated by the force- 
restore method, as derived by Deardorff (1978) and recently modified by Lin 
(1979), i.e.. 


a 


dT 

dt 




P . c -"3- 
sl si 1 


2ir(Tg-T) 

T 1 


( 1 ) 


and 


dT 

dt 


G/(p ^s2 d 2 ) 


( 2 ) 


where a = 1 + >/ ir S/dJ and the ground heat flux 


G = R - XE - H 
n as 


( 3 ) 
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The soil moisture is governed by 


de 3 
J 1 It 


a (I - ~)+ a (X - u.) - q 

o' c t\ * c C T 1 


s s 


w 


12 


for the surface layer or upper root zone and 


do. 


h dt ■ °C U 2 + q 12 " q 23 


(4) 

(5) 


for the intermediate layer or lower root zone, whiles 3 = constant. No 

surface and/or subsurface runoff is allowed in the tests. All rain and 
irrigated water are allowed to infiltrate into the layers. The vegetation 
fraction cr c was set to unity for the vegetation test, while the bare soil 
fraction a s was set to unity for the bare soil tests. The evapotranspi ration 
is parameterized by scaling from the potential evaporation 


E 

P 


- p c V 
a m a 


u*(y - , a ] 


( 6 ) 


as follows. For bare soil 


E 

s 


3 E 
s p 


3 g * min (a 6^,1) 


(7) 


where the normalized available soil moisture is defined as 




e 

pwp 




e ) 

pwp 


i = 1,2 


( 8 ) 


The soil at the Arizona site is Adelanto loam with a = 1.95 and at the 
Georgia site is Cecil sandy loam with a = 4.63, and in both cases, b = 3. 
The a- and b- values were determined by calibration. For a vegetated surface, 
the uptakes of moisture from the upper and lower root zone are parameterized, 
using the concept of soil water stress and the unstressed transpiration. The 
transpiration is the sum of the uptakes 


and 


e t = 


u, 4- u r 


{ 



if e , > 1 

al — 


(9) 


(10a) 


106 


or 


where 


e i¥u- 

If 

1 > 

°al 

E U * 

if 

®a2 

> 1 

Wu • 

if 

1 > 

6 a2 


A . ,0ai/r 

^ = mia [ /0 at 

, 1 

1 



(1 Ob) 


(IT) 


The threshold value of the normalized available soil moisture depends on the 
atmospheric demand such that 


6 at " <CE u /d6 aJ k (12) 

For the Arizona winter wheat, k = 0.8, C = 35 h, and d = 50 cm were used. 
Finally, the unstressed transpiration can be related to the potential 
evaporation in Equation (6) by 


E = aE (13) 

u p 

Where can be estimated from the type of vegetation canopy, e.g. , a = .8 for 
crops. 

The soil moisture flux q-j j in Equations (4) and (5) is the average 
flux in a layer in the neighborhood of the interface, in which the soil 
moisture content of the profile changes from 0i to 9j, i.e. 

q ij = D ij ( 0 i “ 9 j> + K ij » 1 = 1, 2 and 2, 3 (14) 

To evaluate 'Dij and Kij, the soil moisture diffusivity and conductivity are 
assumed as an exponential function of 0, thus 


D±j = (D(0 i ) - DC* Mfdy In { D(G i ) /DC0 j )} 


(15) 


Similarly, 


K ij " {(0 i ) ~ K(0 j )}/ln{K(0 i )/K(O 1 ) } (16) 

It is found that d^j can be estimated by the following scheme: 


di 2 s 1/2 (d] + da) and d23 = 1/ d 2 


(16a) 


if 


Q , > 0 ? > 0 


1 > u 2 > u 3 or 0 i < 0 2 < 0 3 


or 


d] 2 " 1 /4 ( d-j + d 2 ) and dg 3 s 1/4 d 2 


(16b) 


otherwi se. 

The GHM was forced by the observed atmospheric data consisting of 
specific humidity, wind speed, net radiation, and air temperature, and the 
ground temperature and soil moisture were also initialized by the data as shown 
at the upper four graphs in Figs. 2-4. In the Arizona experiments, the plot 
was irrigated prior to the initial time as indicated by the initial moisture 
content in Layer 1. There were two rain events in the Georgia experiment: 
9.7 and 4.1 mm during the 5th and 6th day, which substantially increased the 
soil moisture in Layer 1, as seen in Fig. 3. Calculated results of 0j , ©2» 
Tg, and Ea were compared with the observed data in Figs. 2-4, except that Ea 
was not measured at the Georgia site. All calculated soil moisture contents 
agree very well with the observed data for the weighted-average values in 
each layer. In Figs. 2 and 3, the bare soil ground surface temperature was 
calculated on the basis of <5- 2 and 1 cm, while the observed data were taken 
at a depth of 1 and 0.5 cm for the Arievna and George site, respectively. In 
Fig. 4, the observed canopy temperature was used for comparison. The calcu- 
lated evaporation for the bare soil at the Arizona site in Fig. 2 agrees 
with the observed data when the soil is wet and there is only an overall 
agreement when the soil becomes dry (note the missing data for the 2nd and 
3rd day). In the Georgia experiment in Fig. 3, the evaporation was not 
measured; however, the calculated evaporation was nearly equal to the poten- 
tial evaporation during 2 days after the rain events. The transpiration 
from the winter wheat crop in Fig. 4 shows .a good agreement for the whole 
experiment (note the missing lysimeter data on the 6th day). 

The tests of the GHM using these field-observed data sets have 
provided a useful means to evaluate the parameterizations of soil moisture, 
ground surface temperature, and evapotranspi ration. Good agreement was found 
between calculated and observed values. Further testing is being conducted 
for the parameterization of vegetated surface, using observed data in a forest 
site. 
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NOTATIONS 


a,b,C,k 

C 


m 


D 12' D 23 

d 

E 


E ,E 
p’ s 


G 

H 


I ,1 

s’ c 


si 

K 12’ K 23 

P 

p a 

q*,q a 

q 12’ q 23 

R 1’ R 2 

R 


R ,R 
o’ s 


constants 

moisture transfer coefficient 
soil specific heat 
moisture transfer coefficient 
sum of d] and d 2 

depth for 0],0 2 , diurnal and annual temperature wave 

total evaporative flux, E a * E $ + U] + U 2 

potential and bare soil evaporation 

unstressed transpiration 

ground heat flux 

sensible heat flux 

infiltration on soil and under vegetated surface 

infiltration capacity 

infiltration from snow or ice melt 

soil conductivity 

precipitation or throughfall 

air surface pressure 

specific humidity at saturation and at reference height 
soil moisture flux 

fraction of plant roots 
net radiation 

surface and subsurface runoff 


T depth-average ground temperature 

T a air temperature at reference height 

ice temperature 
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NOTATIONS (Continued) 


g 

U 1» U 2 

V 

a 

e 


^vi 

0 l'V°3 

9 a 

0 fc’°pwp 


p a ,p w' p s 


S* c 


Y 

5 


cell -average ground temperature 

moisture uptake by plants 

wind speed at reference height 

scaling factor for bare soil 

scaling factor for evapotranspiration 

volumetric soil moisture content 

available soil water content 

soil moisture at field capacity and permanent 

wilting point 

latent heat of vaporization 
density of air, water, and soil 
fraction of bare soil and vegetated surface 
diurnal period 

depth correction factor for calculating Tg 
ratio of E u /Ep 

thickness of thin surface layer for calculating Tg 
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INTERANNUAL VARIATIONS OF SEA LEVEL PRESSURE 
(R. Godbole) 


In present-day climate modeling, attention is focused on the use 
of climate models for short-range climate predictions on a time scale of a 
season or a year. An important question to be addressed in this connection 
is the predictability of observed variations in monthly or seasonal means 
from one year to another. A logical prerequisite, therefore, is to examine 
and document these observed year-to-year variations. 

The pu 1 * ,ose of this paper is to examine the spatial patterns of sea 
level pressure and their interannual variations during the past 16 years. 
The period covers from 1961 through 1976 and the region considered is the 
whole globe. January and July monthly mean sea level pressure data for 
individual stations have been plotted and subjectively analyzed. The data 
were extracted from NCAR data sets. 

Figs, la and lb show the 16-year mean sea level pressure for 
January and July, respectively. The major centers of action such as the 
Aleutian low, the Icelandic low, the Siberian high in January, the monsoon 
low over Africa and Asia in July and also the subtropical highs in both 
January and July have the right locations and intensities when compared to 
any corresponding climatological or long-term mean sea level pressure field 
(Shutz and Gates, 1971; Crutcher and Meserve, 1970). 

Fig. 2 represents a typical example of the variation of sea level 
pressure for January of two different years: 1974 and 1976. It is observed 
that in 1974 (Fig. 2a), the Icelandic low is an intense, huge system apparently 
dominating the circulation of the northern hemisphere. The Aleutian low is 
weak and the Siberian high, although strong, is broken into small circulation 
cells. In 1976, on the other hand, the Icelandic low is weak having a second- 
ary center lying over northern Europe. The Siberian high shows a single and 
relatively weak cell whereas the Aleutian low shows intensification. Similar 
interannual variations are noticed for July also but their amplitude is 
considerably less than that observed in January. This is clearly borne out 
from Fig. 3 which shows year-to-year variations of central pressure of a 
major system in January and July. For an Icelandic low in January, the 
maximum range of variation of central pressure is about 32 mb whereas for a 
monsoon low over Asia in July, the maximum range is hardly 5 mb. 

Fig. 4 shows the range of variability in the location of the major 
centers of action for January and July. The results again point out that the 
variability in January is more than the variability for July. Also, these 
year-to-year displacements are stretched more in east-west rather than 
in north-south extent. 

/ Table 1 presents anomalies of sea level pressure for some of the 

major systems in January. The anomalies were obtained as the difference 
between central pressure in an individual year and that of a long-term mean 
(Fig. la) for a given major system. The results reveal a very interesting 
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ytApis 

Year-to-year variation of central pressure (P-1000) for 
Icelandic low in January (solid line) and for monsoon low 
in July (broken line). 








Table 1. Monthly mean sea level pressure anomaly (P-F) for Icelandic low, 
Aleutian low, and Siberian high in January. 


Year 

Icelandic 

low 

- strong 
+ weak 

Aleutian 

low 

- strong 
+ weak 

Siberian 
low 
- weak 
+ strong 

1961 

- 8.0 

- 8.5 

+ 1 .0 

1962 

- 9.0 

+ 0.5 

+ 1.0 

1963 

4-11.0 

- 3.8 

+ 4.8 

1964 

+ 0.0 

- 8.0 

+ 4.0 

1965 

+ 0.0 

- 1.5 

+ 1.0 

1966 

- 2.0 

-10.1 

+ 0.0 

1967 

+ 6.3 

- 4.8 

+ 4.0 


+ 2.2 

+ 2.6 

+ 0.0 


+ 5.1 

- 5.5 

+ 5.1 


-10.1 

+ 2.1 

+ 2.6 


- 6.0 

+ 7.5 

+ 4.0 


-12.8 

+ 3.5 

+ 2.0 


- 8.6 

+ 2.3 

- 1.0 


-20.6 

+ 0.5 

+ 3.0 


- 5.8 

+ 1.1 

■f 1.0 


+ 0.4 

- 7.8 

- 1.2 
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feature. It is seen that, with the exception of 2 to 3 years, the Icelandic 
low bears an inverse relationship with the Aleutian low. This means that, 
in a given year, if the Icelandic low is strong, the Aleutian low is weak 
and vice versa. Further, the behavior of these two lows seems to be somewhat 
associated with that of the Siberian high. The relationship, however, is 
not as strong as that noticed between the two lows. Nevertheless, one cannot 
fail to notice that in a majority of cases, a weak Icelandic low tends to 
associate with a weak Siberian high whereas a weak Aleutian low tends to 
associate with a strong Siberian high. A somewhat similar relationship 

appears to exist between the subtropical highs also. The findings have a very 
important implication: these seasonal centers of action which act as internal 

forcing in a climate model dynamics have an impact not only restricted locally 
where they are produced but rather their impact is global in nature. 
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PRELIMINARY RESULTS OF A JANUARY SIMULATION 
WITH AN IMPROVED VERSION OF THE GLAS MODEL 

(J. Shukla, Y. Sud, and E. Sabatino) 


Several aspects of the weaknesses and strengths of winter and summer 
simulations with the GLAS model are described in a recent paper by Halem et al . 
(1979), which points out that one of the weaknesses of the model is its lower 
evaporation over the oceans. This, along with smoothing of orography over 
high mountain peaks, was suspected to be responsible for weaker amplitudes 
of stationary waves in these simulations. We have introduced several modifi- 
cations in the earlier version of the GLAS model. These were comprised of 
modified boundary layer calculations which led to an improvement in the fields 
of boundary layer fluxes and a stable calculation of surface temperature. 
This parameterization is discussed in a paper by Sud and Abeles (1979). In 
addition, the ground albedo over snow-ice covered surfaces of the earth are 
prescribed from an analysis of satellite observation of Matson (1978) and 
Rayner and Howarth (1977). This snow-ice data is combined with the bare 
land albedos of Posey and Clapp (1964). Furthermore, the model has been 
initialized for a long climate run by making the permanent ice depth equal to 
3 meters and transient snow depth of 0-2 cm, depending upon its reflectivity. 
Sea surface temperatures from a 2° x 2° GFDL data set were interpolated to 
the 4° x 5° grid of the GLAS model. All the prescribed boundary conditions 
are now interpolated on a daily basis. 

Furthermore, as a result of detailed examination of the code and 
the diagnostics of various calculations, a few additional modifications have 
been made. These consist of the following: elimination of Shapiro-filter on 
humidity field, maintaining conservations in strapping and unstrapping of 
cloud layer in the model, calling short wave radiation every half hour while 
maintaining long wave calls at 5-h frequency but with net radiative tempera- 
ture modification to be made every half hour using the saved fluxes of long 
wave radiation, changes for correction of radiative heating rates over the 
mountains, and correcting a coding error in the interpolation of level -8 
winds in the split-grid regions. 

Preliminary results of first 31 days which consist of simulated 
month of January (starting from the initial conditions of 1 January 1975 
from NMC analysis) are compared with the simulation with the earlier model 
run of Halem et al . (1979) of which February simulation has been presented 
in their paper. Some basic fields which show an impact from these changes 
are presented and compared with the earlier run. 

Analysis of Results 

An important result from these modifications was the improvement of 
one of the most serious deficiencies in the model, i.e., the lower amplitude 
of the stationary component of its general circulation. It is well known that 
the stationary waves are forced by diabatic heating sources and orography. 
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Since the earlier model produced an evaporation field which was significantly 
less than the observed over the oceans {Charney et al . , 1977)* an improvement 
in diabatic heating through increased boundary layer fluxes of evaporation 
and sensible heating over the oceans was expected to improve the stationary 
component of general circulation. This has happened in the new model 
simulation. 

The following results compare the modified model with its earlier 
version. Figure captions (a) and (b) are used for the modified and earlier 
model results, respectively. Figs, la and lb show the sea level pressure 
maps. Fig. la has better subtropical highs and a pronounced Siberian high. 
Also, the low over the coast of Japan is much better defined and the sea level 
pressure fields off the coast of Greenland have a better structure. 

Figs. 2a and 2b show zonal plots of evaporation over the land and 
oceans in the two simulations. Increased evaporation over the oceans and 
lower evaporation over land in the modified model run are seen. This corrects 
a well known weakness of the model. Figs. 3 and 4 show the zonal ly averaged 
evaporation and precipitation fields. There is an increase in evaporation 
and precipitation in the new model. 

An analysis of amplitude of stationary component of circulation 
geopotential height fields are presented. Tables la and lb show the variance 
of geopotential height fields at the mandatory pressure levels. The values 
of variances for the new run compare well with the observations (as shown 
in Table 1c) of Oort and Rasmusson (1971) and show a considerable improvement 
over those of the old model shown in Table lb. In Fig. 5 the plots of variance 
of geopotential height fields at 500 mb are shown. The new model results are 
in good agreement with observations. The peaks of variances in the baroelinic 
regions improve by a factor of two over the old model. 
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* Table 1c. Variance of geopotential height resulting from stationary eddies [Z* 2 ] for the month. 

(units in 10 2 gpm 2 ) 
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January 




(10 2 gpm?) 



P(mb) 

10S 

5S 

EQ 

5N 

10N 
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20N 

25N 

30N 

35N 

40N 

45 N 

50N 

55N 

60N 

65N 

70N 

75N 

50 

6.4 

6.4 

13.6 

18.7 

22.2 

25.8 

25.6 

29.2 

32.3 

41.0 

54.1 

85.9 

111.4 

150.4 

178.3 

189.8 

145.0 

71.8 

100 

1.9 

3.3 

6.0 

7.3 

9.0 

11.7 

13.5 

14.6 

21.1 

41.6 

75.8 

136.2 

191.0 

223.7 

257.4 

253.1 

209.8 

124.7 

200 

1.6 

1.6 

2.9 

4.3 

7.0 

11.7 

16.9 

20.0 

31.4 

71.7 

135.8 

219.3 

258.4 

243.6 

206.5 

164.1 

122.1 

72.5 
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.2 
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2.3 

4.0 

6.3 
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14.3 

36.6 
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168.2 
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.5 
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1.2 

1.6 

2.4 

3.8 

5.8 

9.3 

20.3 

45.9 

81.3 

122.1 

137.9 

122.0 

98.2 

79.5 

64.5 

42.7 


* As taken from Table Bllb of Oort and Rasmusson (1971). 
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Fig. 3. Zonally averaged evaporation for the modified and earlier GLAS models. 
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Variance of 500 mb geopotential height for the 
modified and earlier GLAS models and observations. 


MOISTURE AND ENERGY BALANCE STUDY FOR A SET OF WINTER 
AND SUMMER SIMULATIONS WITH THE GLAS MODEL 

(Y. Sud) 


A heat budget analysis of the earth atmosphere system yields heat 
balance fields separately at the surface of the earth, within the atmosphere, 
and at the top of the atmosphere. Basically the analysis is helpful in under- 
standing the nature of energy and moisture storage in the atmosphere, trans- 
fers at the boundaries of the atmosphere, and transport within the atmosphere, 
'i he analysis is useful in correlating the dynamical processes with the energy 
transfers and transformations. In a GCM simulation, it has yet another very 
useful purpose, i.e., checking the conservation of energy and moisture globally. 
In the GLAS climate model, lack of conservation in the Shapiro filter is well 
known. Certain other ad hoc adjustments such as zeroing out of negative humidity 
and making the minimum specific humidity equal to 10" 8 also contribute to 
lack of conservation. There are similar nonconservative simplifications in 
the physical calculations. Whereas in the next climate model all these incon- 
sistencies will be eliminated to get even better conservation, the present 
model does yield some unaccounted energy and moisture in the atmosphere. 

The circulation of global atmosphere is supported by thermal forcing. 
There are two basic sources of thermal forcing on the time scale of a season, 
which is 3 months or less in these studies: solar insolation and the fluxes 

of heat and moisture from the surface of the earth including oceans. Oceans 
are represented by a prescribed water surface temperature which varies steadily 
over a month. This implies that the oceans are infinite sources of sensible 
and latent energy. In a realistic interaction between the oceans and the 
atmosphere, a reasonable balance between the net radiation into the oceans and 
the fluxes of sensible and latent energy from the oceans is vital. On the 
surface of land, a simple conservative slab model is used. The change in soil 
temperature and soil moisture is determined by the net heat and moisture flux 
(precipitation-evaporation-runoff) to the soil. Over snow and ice on land, a 
part of the heat flux to the land may be used in melting ice or snow. This 
flux is lost for the surface heat balance through runoff. All energy and 
moisture balances are based on the conservation principle and may be derived 
from energy and moisture-conservation equations. A detailed discussion of 
the GLAS global circulation model may be found in Somerville et al . (1974) and 
Halem et al . (1979). 

The GLAS model climate simulations generate diagnostics of energy 
and moisture fluxes at the top and bottom boundaries of the atmosphere. In 
conjunction with this, the initial and final states of the atmosphere for the 
simulation period may be used to obtain a complete energy and moisture cycle 
analysis of the earth-atmosphere system. The net surface balance can also be 
used to give the energy flux to the ocean. From all this data, calculations 
can be made to obtain energy transport in the atmosphere as well as the ocean. 
Similar calculations for moisture field yield moisture transport in the atmos- 
phere. This study examines these quantities in great detail. Presented here 
is an analysis of a set of three summer and three winter runs. Each set 


consists of a control run from the observed initial state of the atmosphere 
plus two predictability runs in which only the wind fields at the initial 
time are perturbed by 3 m sec" 1 . The heat balance is obtained in the form 
of two-dimensional, zonal, hemispheric and global distributions of the sought 
quantities. Also, energy balance at the top of the atmosphere, energy balance 
at the surface of the earth, and moisture and heat balance in the atmosphere 
are obtained. Taking into account the storage of moisture and other forms 
of energy, the transport of heat and moisture in the atmosphere and transport 
of heat in the oceans are calculated. 

Governing Equations 

The energy and moisture flux diagnostics for the atmosphere are 
generated for an overall vertical column of air. There is no information 
about breakdown of these fluxes for individual a-levels of the model. Water 
balance may be calculated from the balance between the fluxes of evaporation 
and precipitation. The following fluxes are written out as diagnostics every 
6 h and their values are accumulated means over various model time steps in 
a given time interval. The units of each quantity are cal cm -2 day - '. 

F s = Upward evaporation flux at the surface 

P s * Precipitation flux to the surface 

H s = Upward sensible heat flux at the surface 

Sy = Incoming short wave radiation at the top of the atmosphere 
S a = Absorbed short wave radiation by the atmosphere 
S s = Upward long wave radiation at the surface 

Ly » Upward long wave radiation at the top 

L s = Upward long wave radiation at the surface 

Hp = Heat storage in melting/freezing water and heat conduction through 

sea ice. 

The above diagnostics along with the initial and final states of 
the atmosphere are used to obtain a complete heat and moisture balance. The 
basic equations follow: 

The total change of a conservative variable, s, representing energy 
or mass as the case may be, can be written as, 

( e . ‘ E .)/At 
at f i 


A1 so, 

de _ 9£ + V. (evl - JL (ecn , 
dt at ' ?o 

ae a a 

substituting for from (1) and taking^ = -9JL , one obtains 
at dt dt 


( 1 ) 

( 2 ) 


ds = ( e - e )/ At + V . (ey) 
dt f r 
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Obviously = 0 for a column of air. Here ds/dt is the source 
term. It represents net flux from the source, s, which may be calculated from 
the diagnostics as follows: 


For heat flux absorbed by the atmosphere, 

is * S + H + P + (L - L ) , 

dt ass s T 

For moisture flux to the atmosphere, 

ds = F - P . 
dt s s 

For energy flux to the oceans/land, 


4i s s - L - F -H + H 
dt s s s s p 


( 4 ) 

( 5 ) 

( 6 ) 


Similarly, e f or £ ] can be calculated at the final and initial time by: 
vertical integrals of enthalpy plus kinetic energy of the air, and vertical 
integrals of specific humidity for moisture in the air. The time integration 
of product C z Tg summed over land gives the heat absorbed by the ground. Here 
C z is the thermal capacity of the ground and T g is the ground temperature. 

Global Integrals 


The adyection calculation in Eq. (3) applies to a finite area. If 
the equation is integrated over the entire globe, the advection integral must 
vanish. Therefore, 


-ft J g|- dA = jj- J [(e f -e i )/A t]dA ' ^ 

A A 

Since both the RHS and LHS of Eq. (7) can be independently calculated over the 
globe, the two sides should be equal in a conservative model. Any difference 
between the two represents lack of conservation. This difference is obtained 
for the energy and moisture balance in the atmosphere as well as the net heat 
flux to the oceans. Satellite observations analyzed by Barrett (1974), and 
Vonder Haar and Suomi (1971) show that global integrals of energy flux to the 
earth-atmosphere system should be vanishingly small. This implies thermal 
equilibrium of the oceans with the atmosphere which is a consequence of oceans 
attaining equilibrium with the atmosphere for millions of years. 

Analysis of Results 

From the monthly average of flux fields, we obtain their zonal, 
hemispheric and globally averaged distributions. Tables la and lb show the 
hemispheric and global heat balances for February and July simulations. The 
mean fields and their standard deviations between the three runs for each 
hemisphere as well as over the entire globe are also shown. From the fluxes 
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Table la. February 1976 energy balance (cal/cm 2 /day). 


Analysis of heat budget fields 


S. hemisphere N. hemisphere 



Top boundary of atmosphere 

Incoming short wave flux 
Outgoing short wave flux 
Outgoing long wave flux 


Within the atmosphere 

S.W. absorbed 
L.H. release by 
precipitation 

S.H. increase (atmosphere) 
KE increase (atmosphere) 

LE increase (moisture) 


S.W. absorbed 
L.W. outgoing 
S.H. fl ux outgoi ng 
L.H. flux outgoing 


Net e 


La nd 
Ocean 

All surfaces 

At the top (atmosphere) 


Unaccounted quantity 

Heat in the atmosphere 
Moisture in the atmosphere 


Mean SO 


915.90 0.00 
284.50 1.35 
429.70 2.10 


125.60 

164.20 

-9.40 

-0.15 

-0.80 


505.70 1.78 
113.80 0.82 
48.00 0.53 
145.40 0.67 


198.40 1.00 
-201.60 0.56 


Mean SD 


565.80 0.00 
170.60 1.03 
402.50 1.80 


78.40 0.14 
134.20 1.00 

3.75 0.80 
- 0.01 0.11 
1.11 0.80 


316.70 0.66 
115.90 0.39 
70.00 0.67 
145.00 0.51 


-14.08 0.84 
8.34 0.18 



Global 

Mean 

SD 

740.90 

0.00 

227.60 

1.19 

416.60 

1.95 

102.00 

0.33 

149.20 

1.10 

-2.84 

0.57 

0.07 

0.08 

0.16 

0.22 

411.20 

1.22 

114.90 

0.60 

59.00 

0.60 

145.00 

0.59 

10.76 


136.30 

0.32 

92.18 

0.92 

-96.65 

1.37 

11.23 

0.32 

-4.14 

0.12 


0.31 0.15 


0.30 0.18 


0.30 0.16 


i 




















Table lb. July 1974 energy balance (cal/cm 2 /day). 



Analysis of heat budget fields 

Items 

S. hemisphere 

N. hemisphere 

Global 


Mean 

SO 

Mean 

SD 

Mean 

SO 

Top boundary of atmosphere 







Incoming short wave flux 

442.30 

- 

948.90 


695.60 


Outgoing short wave flux 

129.50 

1.37 

307.49 

2.95 

218.52 

2.06 

Outgoing long wave flux 

388.70 

0.61 

438.26 

1.06 

413.46 

0.84 

Within the atmosphere 







S.W. absorbed 

66.00 

0.08 

142.79 

0.23 

104.40 

0.16 

L.H. release by 







precipitation 

95.10 

1.30 

135.90 

1.75 

115.50 

1.53 

S.H. increase (atmosphere) 

-7.89 

0.70 

18.18 

0.19 

5.14 

0.38 

KE increase (atmosphere) 

0.62 

0.06 

0.02 

0.01 

0.32 

0.03 

LE increase (moisture) 

-3.28 

0.62 

7.11 

0.09 

1.92 

0.46 

Bottom boundary of atmosphere 







S.W. absorbed 

246.75 

1.29 

498.61 

2.83 

372.68 

2.06 

L.W. outgoing 

148.44 

0.67 

178.25 

1.10 

163.35 

0.89 

S.H. flux outgoing 

41.55 

1 .05 

43.08 

0.91 

42.31 

0.98 

L.H. flux outgoing 

112.03 

1.14 

112.94 

0.79 

112.48 

0.97 

Net energy absorbed 







La nd 



_ 


5.63 


Ocean 

- 

- 

- 

- 

75.27 

4.06 

All surfaces 

-55.27 

2.41 

164.35 

0.84 

54.54 

1.62 

At the top (atmosphere) 

-75.90 

0.97 

203.15 

1.85 

63.62 

1.47 

Unaccounted quantity 



• 




Heat in the atmosphere 

- 

- 


_ 

14.68 

4.28 

Moisture in the atmosphere 

- 

- 


*• 

-13.70 

3.96 

Planetary albedo 

0.29 

0.003 

0.32 

0.003 

0.31 

0.003 
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at the top and bottom of the atmosphere, the net energy balance at the top and 

bottom of the atmosphere is calculated. The surface energy balance is further 

subdivided into its land and ocean distributions. Land is treated as a slab 
of finite but small heat and moisture capacity. Therefore, averaged over a 
month, the surface heat and moisture balance over land should be negligible. 

On the other hand, over the oceans the net surface fluxes may be substantial 

depending upon the region. However, in a realistic model, the net globally 
averaged surface flux over the ocean must be small. This is found in observa- 
tions (Vonder Haar and Suomi, 1971). This result also must follow as a 
natural consequence of oceans equilibrating with all its environments over 
millions of years. In fact, this is one of the basic assumptions in simpler 
climate models which neglect oceanic storage (e.g., Sellers, 1967). However, 
this is pot the case with our model simulations. 

There is a large net globally averaged heat flux sinking into the 
oceans both in the February and duly simulations. This situation in turn 
leads to reduced outgoing long wave flux at the top of the atmosphere. This 
reduced flux must be supported by lower upper level temperatures of the 
simulation model which is indeed one of the weaknesses of our simulation 
model. Thus, too cold upper atmosphere simulated by the model seems to be 
consistent with sinking of heat in the oceans. However, it is not clear how 
the two effects are dynamically related. We know that the present climate 
model oceans do not evaporate sufficiently. Past attempts to increase 

evaporation by an ad hoc increase in the boundary layer transport coefficients 
in the model were not successful because the final equilibrium of ocean- 
atmospheric interactions brought the two to a state which leads to very 
little change in the fluxes of sensible heat and evaporation over the ocean. 
All these findings point out that the main weakness may be in the upward 
transport of energy. In whatever way this discrepancy is resolved, it 
appears that cooling of upper model levels and sinking of net energy in the 
oceans are manifestations of the same problem. 

In Table 2 computed upwelling long wave fluxes are shown for two 
cloud-free columns of air which have temperature and humidity profiles for 
the mean monthly February zonal ly averaged at SON for the model simulations 
and analyzed observations (Oort and Rasmusson, 1971). Notice the model 
temperature profile has typically colder upper atmosphere and yields 8.2 per- 
cent lower outgoing long wave flux at the top of the atmosphere. The re- 
duction in cooling is of the order of heat sunk in the oceans. The heat 

balance is shown schematically in Figs. 1(a) and 1(b) for the data in 
Tables la and lb for February and July, respectively. 

Concl usions 

The heat balance analysis of February and July simulations with the 
GLAS model reveals that the GLAS model energy and moisture balance fields are 
in good agreement with observations and present a consistent picture. Con- 
servation is reasonable but is not within the limits of round-off errors. 


Table 2. Computed upwelling long wave fluxes for two cloud-free 
columns of air for February, zonally averaged at BON. 



Model simulation 


Observed - * - 

Pressure 

(mb) 



Temperature 

(°K) 

Humidity 

(9/kg) 

Upward long 
wave flux 
(cal/cnf*/day) 

Temperature 

HO 

Humidity 

(g/kg) 

Upward long 
wave flux 
(cal/cm 2 /day) 

65 

203.6 

0.02 

415.8 

218.2 

0.02 

449.9 

175 

205.9 

0.12 

404.2 

217.2 

0.15 

427.3 

285 

214.3 

0.31 

391.5 

220.2 

0.35 

406.8 

395 

228.7 

1.34 

372.8 

232.2 

2.50 

383.3 

505 

242.0 

4.10 

339.1 

244.2 

4.00 

343.2 

615 

252.3 

7.90 

300.6 

253.2 

8.00 

309.4 

725 

259.5 

11.90 

264.8 

259.2 

12.50 

274.3 

835 

264.5 

16.80 

233.5 

264.2 

17.50 

244.7 

945 

267.0 

26.40 

203.6 

268.2 

27.80 

216.8 

1000 

268.3 

30.00 

175.2 

270.2 

30.00 

186.9 


' Observed data represents winter season mean taken from Oort and Rasmusson (1971) 
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The model loses .02 mb/day of air mass which is within the limits of round-off. 
However, making up this loss by adding mass uniformly itself contributes to 
lack of conservation. The major factors for the lack of model conservation 
are: use of nonconservative Shapiro filter, model round-off, and some other 
nonconservative calculations discussed in the introduction. Another major 
weakness of the model is its sinking of large amounts of energy in the oceans. 
The next climate model which uses a modified boundary layer to enhance fluxes 
of moisture and sensible heat over the oceans (Sud and Abeles, 1980) is ex- 
pected to correct this situation. The study reveals that upper level cooling 
of air may be a consequence of lack of surface energy balance in the model 
brought about by adjustment through dynamical transport processes. If it 
turns out that this is indeed the case, then the finding Is important because 
It points out that the boundary layer processes can affect the upper levels 
of the atmosphere on a climate scale. However, this result is merely an 
indirect inference and needs further evidence. 
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ARCTIC AND ANTARCTIC CLIMATOLOGY OF THE GLAS GCM 
(G. F. Herman and W. T. Johnson) 


A detailed study of the high latitude climatology of the GLAS GCM was 
conducted in order to assess the reliability of the model for high latitude 
climate studies, such as those involving changes in atmospheric carbon dioxide 
or sea ice boundaries. 

The model -generated January climatology was constructed from eight 
simulations (Runs No. D058, D077, D084, D091 , D092, D093, D094, D095), while 
the July climatology was based on seven simulations (Runs No. D115, D116, D119, 
D120, D126, D127, D129). These runs differed from one another by either small, 
randomly distributed values in the initial conditions, by variations in sea 
surface temperature or sea ice margin, or by minor variations in model physics. 
The selection of runs involving both different initial conditions and different 
boundary conditions is an attempt to represent some of the processes in nature 
that cause one January (or July) to differ from another. 

The January climatology of the model was computed by forming the 
average of eight 30-day simulations, which in turn were each formed by averaging 
the model -generated fields sampled at 12-h intervals. A comparison is made 
between these results and the NCAR data set made available by Jenne (1975). 

Fig. 1 illustrates the distribution of sea level pressure for January 
in the northern hemisphere as obtained in the model and given by the NCAR 
climatology. The dominant features of the observed circulation are the major 
low pressure systems over the North Atlantic and North Pacific, the broad and 
intense Siberian anticyclone, and a somewhat weaker anticyclone over northwest 
Canada. The Icelandic low in the model's climatology is in excellent agreement 
with observations, but appears to be about 4 mb too deep. The model's Aleutian 
low is not quite as extensive as observed, and its center is located about 
15° longitude too far to the east. Ire major deficiency of, the model's high 
latitude climatology is its failure to represent in any realistic way the 
Siberian anticyclone. In nature, this high pressure center extends over a vast 
region of the central U.S.S.R. from the Caspian Sea to the Sea of Japan, with a 
central pressure exceeding 1032 mb. The high in the model is confined to 
northwest Siberia, with,central pressures on the order of 1024 mb. 

The distribution of pressure during the northern hemisphere summer 
is shown in Fig. 2. During the summer, the observations show that the winter- 
time lows over the oceans have been replaced by broad regions of high pressure, 
while weak thermal lows dominate the continents. An area of weak low pressure 
is found in the Central Arctic. In the model climatology, the maritime high 
pressure areas tend to be too weak by about 4-6 mb, and the continental lows 
too deep, particularly over Asia. A weak ridge of high pressure spans the 
Central Arctic. 

Sea level pressure over Antarctica is not well defined because of 
the difficulty both in the model and in the observations of reducing surface 
pressure to sea level. However, a comparison can be made over the peripheral 
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sub-Antarctic oceans. Observations indicate that during the winter (July) 
(Fig. 3), the Antarctic continent is encircled by deep low pressure regions 
found in the Ross, Weddell, and Bel li nghausen Seas. Broad regions of high 
pressure are found in the central regions of the Atlantic, Pacific, and 
Indian Oceans. The July model climatology misses most of these features. 
While the ridge of high pressure in the South Pacific is well simulated, 
that in the Indian Ocean is missing, while that in the South Atlantic is 
underestimated. The stationary low pressure centers around Antarctica do 
not occur in the model, and this points to a problem that the model has in 
generating baroclinic eddies in the high latitudes of the southern hemisphere. 
During the summer (January) (Fig. 4), the high pressure regions over the 
oceans are even more expansive, and a number of low pressure centers continue 
to ring the Antarctic continent. The model does a reasonable job of repro- 
ducing the oceanic highs in January, but is capable of producing only a weak 
low pressure center (996 mb) in the high latitudes of the South Pacific. 

The mean flow at 500 mb during January is shown in Fig. 5. The 
dominant features of the observed January climatology are deep troughs located 
over the Hudson Bay and the Sea of Okhotsk, and a weaker trough over the western 
U.S.S.R. Maximum zonal flow occurs in midlatitudes off the east coast of Japan 
and over the eastern U.S. In the model both the midlatitude jet and the 
Canadian trough are well represented, but the axis of the trough is located 
about 15° longitude too far to the east. A Eurasian trough appears, and also 
is slightly displaced to the east. The most glaring deficiency of the model 
appears to be its failure to reproduce the Pacific trough, which is a diffi- 
culty likely to be associated with problems in treating adequately the Hima- 
layan topography. 

It appears that, at the surface in the northern hemisphere, the model 
does its best at reproducing circulation features over the oceans. This is not 
surprising since the oceans exert such a strong influence on surface features 
and have their temperatures specified in the model. The agreement with obser- 
vations is not as good over land, and it is speculated that this is a conse- 
quence of problems in treating the interaction between the surface radiation 
budget and cloudiness. 

The long wave structure of the model's climatology appears reasonable 
in the western hemisphere but not in the eastern hemisphere. 
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SEA ICE SIMULATIONS BASED ON GLAS GCM OUTPUT 
(C, L. Parkinson, G. F. Herman, M. R. Good, and W. T. Johnson) 


The sea ice cover of the north and south polar seas plays an 
important role in determining the energy balance of the surface through its 
effect on surface albedo and its influence on heat and moisture transport at 
the air-sea ice interface. At the present time, the sea ice coverage of the 
GLAS GCM is a specified parameter, i . e, , it varies according to a predetermined 
climatology and does not alter the value of model -generated variables. 

As a preliminary step in the development of a fully coupled ice- 
atmosphere-ocean model , an "offline" sea ice simulation is being conducted in 
which the atmospheric fields calculated with the GLAS GCM have been used as 
forcing for a sea ice model, but in turn are not affected by the resulting 
distribution of sea ice. 

The sea ice model employed is essentially the thermodynamic portion 
of the Parkinson and Washington (1979) model. It includes a snow layer, an 
ice layer, and a lead percentage within each grid square, The calculation 
involves energy balances at the water-ice, ice-snow, snow-air, and water-air 
interfaces. The incoming short and long wave radiation and precipitation 
amounts are based directly on the GCM output, while the sensible and latent 
heat fluxes are calculated through bulk aerodynamic formulae which use GCM- 
generated air temperatures, specific humidities, and wind speeds. 

The GCM calculation was a 3-month run (GMSF reference No. D1 22) 
based on 00Z 1 January 1975 NMC initial conditions. Initial sea ice concen- 
trations for the same date were estimated from maps of average brightness 
temperature recorded by the Electrically Scanning Microwave Radiometer (ESMR) 
on NIMBUS 5 over the 3-day period, 30 December 1974-1 January 1975. Initial 
ice thicknesses were estimated based on climatology. The initial ice extents 
and thicknesses are shown in Fig. 1. 

The results of the simulation are shown in Figs. 2 through 4. In 
the Arctic, there is some adjustment to the relatively thick initial ice cover, 
and ice in the coastal region melts during January. However, in February and 
March the ice thickens again, as would be expected for that time of year. The 
northern hemisphere ice extents move equatorward from January to March, with 
retreat of the ice beginning by 1 April. 

In the southern hemisphere (Figs. 1-2), the ice melts during January 
and February ( i . e. , during the Antarctic summer). Although the ice thickens 
somewhat from mid-February to mid-March (Figs. 2-3), the ice edge retreats, 
especially in the Bellingshausen and Amundsen Seas. The continued ice retreat 
into March does not correspond well with the observed seasonal cycle as inter- 
preted from the ESMR imagery for 1975. The simulated ice expands outward 
during the latter half of March (Figs, 3-4), and it appears that the simula- 
tion has simply extended the melt season by 1/2-1 month. A likely cause of 
the extended melt season in the Antarctic may be the unrealistically high 
surface temperatures which characterize the model's Antarctic climatology 
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Fig. 1. Initial conditions — ice thicknesses, 1 January 1975. 













Actual thickness and concentration values are listed in Table 1 
for the Weddell Sea, East Greenland Sea, and Central Arctic. The values 
confirm the extension of the melt season into mid-March in the Weddell Sea, 
but also point out the increase in concentration from the start to the 
middle of the month. The full (99.5-100 percent) ice cover in the Central 
Arctic throughout this period is realistic, though thicknesses of 2. 7-2. 8 m 
are about 0,5 m below those typically cited. 


Table 1. Areally averaged ice thickness and concentration. 



Weddell Sea 

E. Greenland Sea 

Central Arctic 

Thickness 

'TMetersJ 




1 January 

1.81 

1.66 

2.73 

15 January 

1.59 

1.51 

2.73 

1 February 

1.35 

1.39 

2.72 

15 February 

1.19 

1.29 

2.72 

1 March 

1.09 

1.23 

2.72 

15 March 

0.99 

1.08 

2.74 

1 April 

1.01 

0.90 

2.76 

Concentration 

(Percent) 




1 January 

60 

62 

5 JO 

15 January 

53 

70 

100 

1 February 

48 

73 

100 

15 February 

45 

74 

100 

1 March 

43 

74 

100 

15 March 

47 

79 

100 

1 April 

67 

79 

100 
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THE EFFECT OF CLOUDS ON THE EARTH'S 
SOLAR AND INFRARED RADIATION BUDGETS 

(G. F. Herman, M-L. C. Wu, and W. T. Johnson) 


A series of numerical experiments were conducted with the GLAS GCM 
to demonstrate the effect that the global distribution of cloudiness has on 
the solar and infrared components of the radiation balance, and in particular 
to determine the relative role of the albedo and greenhouse effects of clouds 
when they are considered in the global aggregate. The complete analysis of 
this experiment has been published elsewhere (Herman et al . , 1980), and some 
aspects of the discussion of the clouds' effect on the earth's radiation 
balance were presented at the Fourth NASA Weather and Climate Program Science 
Review ( cf . Herman et al . , 1979) and will not be repeated here. The focus 
here is on cloud effects on surface and atmospheric temperature, and the 
resultant influence ("feedback") on cloud formation processes. 

In the control run for this experiment, solar and thermal radiation 
calculations were fully consistent and interactive with the model's predicted 
cloud fields. Solar radiation was calculated according to the method of Lacis 
and Hansen (1974), which explicitly treats cloud absorption and multiple scat- 
tering processes. Long wave radiation was computed and summed over ten spec- 
tral intervals using a technique developed by M-L. C. Wu and L. D. Kaplan. 
Cloud solar optical properties were fixed according to cloud type as described 
by Somerville et al . (1974, Table 2), and unit cloud emissivity was assumed 
in the infrared spectrum. Since the model does not treat fractional cloudi- 
ness that occurs on the subgrid scale, clouds were assumed to cover an entire 
4° latitude by 5° longitude grid area. 

In one experiment ("inactive solar"), all cloud formation processes 
and radiative transfer were calculated as they were in the control, except 
that cloud coverage was fixed at zero in the solar radiation calculation. 
Thus, clouds formed, liberated latent heat, and participated in vertical 
mixing processes, but were otherwise transparent ot the streams of solar 
radiation. The second experiment ("inactive thermal") specified zero cloud 
coverage in the long wave radiation calculation, but explicitly treated all 
other cloud processes, including solar radiation. 

The initial state was defined from the 0000 GMT 1 January 1975 
National Meteorological Center (NMC) global analysis, and the integration 
period was 30 days. Conclusions therefore refer to an average January circu- 
lation for a model atmosphere in which the ocean surface temperature, surface 
albedo and sea ice extent are prescribed and non-interactive. Thus, there 
is no means to determine from these experiments whether the annual average 
planetary mean surface temperature would increase or decrease in association 
with global cloud coverage changes. The differences between the control and 
experimental ("transparent" ) simulations represent only the responses of the 
atmospheric and land components of the model's earth-atmosphere system. 
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Analysis of these experiments determined that the radiation balance 
of the January control was more positive (less negative) than the observations 
of Raschke et al. (1973) in the northern and southern hemispheres, and globally 
by 14, 28, and 22 W nr 2 , respectively. This discrepancy is a 1 ikely consequence 
of the model's treatment of thermal absorption in high clouds, or the absence of 
partially clouded areas. When clouds do not interact with solar radiation, 
substantially more energy is available to the system because of the diminished 
planetary albedo, and the balance became more positive (and less negative) than 
the control by 26, 87, and 56 W nr 2 in the northern "and southern hemispheres 
and globally. The balance became more negative (and less positive) in these 
same regions when clouds are transparent to thermal radiation by 35, 45, and 
40 W nr 2 , respectively. 

Under certain circumstances, it is possible to compute the net effect 
of the clouds on the radiation balance simply by combining the radiation loss 
due to the albedo effect, and the gain due to the greenhouse effect. Thus, 
the net cloud effect on the global net radiation balance is -16 W nr 2 , i.e., 
the al bedo effect makes the balance more negative by 56 W nr 2 (39 minus 95), 
while the greenhouse effect makes it more positive by only 40 W nr 2 (39 minus 
-1). Similarly, in the southern hemisphere, a cloud effect of -42 W nr 2 is 
obtained and in the northern hemisphere, +9 W nr 2 . Hence, during the month 
of January, the albedo effect dominates the greenhouse effect in the radiation 
budgets of the globe and of the southern hemisphere. The greenhouse effect 
dominates the budget of the northern hemisphere winter because there is less 
solar radiation available than in the southern hemisphere due to geometrical 
factors. 


It is possible to deduce the net effect of clouds on the net radiation 
as the difference between the albedo and greenhouse effects only if the changes 
caused by the two processes are independent. It is possible to imagine a 
variety of coupling mechanisms or feedbacks that might relate the cloud-solar 
and cloud-long wave processes. For example, ground temperature differences 
induced by solar radiation changes might change the long wave loss at the top 
of the atmosphere through cloud formation processes related to stability and 
cumulus convection. 

Table 1 illustrates the effects of cloud-radiative interactions on 
surface and atmospheric temperatures. Values are spatial averages with appro- 
priate areal weighting. The surface temperatures over oceans do not change 
because ocean tempertures are specified in the model. When clouds are not 

present in the solar radiation calculation, surface temperatures over land 
(Row 1) increase in the northern and southern hemispheres and globally by 2, 
3, and 3°C, respectively, due to the increased direct solar heating of the 
ground. When clouds are absent in the long wave calculation, back radiation 
from the atmosphere decreases substantially and hence averaged surface tem- 
peratures decrease by 5, 2, and 3°C. Spatially averaged atmospheric tempera- 
tures behave in the same fashion, over both land and ocean; i.e., temperatures 
increase by 1-2°C without solar interactions and decrease by 1-2°C without 
thermal interactions. 


Table 1. Spatially averaged temperatures for radiative-interactive 
and non-i nteractive cloudiness (°K). 




(A) 


(B) 


(C) 


January 

Inactive thermal 

Inactive solar 


control 

mi nus 

control 

minus 

control 


Land 

Oceans 

Land 

Oceans 

Land 

Oceans 

Surface temperature 







N. hemisphere 

271 

293 

-5 

0 

2 

0 

S. hemisphere 

285 

291 

-2 

0 

3 

0 

Global 

275 

292 

-3 

0 

3 

0 

Average atmospheric 
temperature 



• 




N. hemisphere 

240 

252 

-1 

-1 

2 

1 

S. hemisphere 

249 

253 

-1 

-1 

2 

2 

Global 

243 

253 

-1 

-2 

2 

1 


The degree of coupling between the radiation and the cloudiness is 
illustrated in Table 2, which shows the differences between the cloudiness 
calculated in the control and in the transparent cloud simulations. As a 
measure of “cloudiness", cloud frequency is used, which is defined as the 
fraction of the total integration time during which a specified cloud type 
occurred at a grid point. 

The supersaturation (or stratus) cloud frequencies increase by 
1-2 percent when cloud-long wave interactions are absent because of the 
cooler atmospheric temperatures ( cf . Table 1). The cloud frequencies de- 
crease by comparable amounts due to warmer atmospheric temperatures when 
solar interactions are absent, except over continental regions of the 
southern hemisphere. Evidently, there is a negative correlation in the 
model between the occurrence of layered cloudiness and average atmospheric 
temperatures. 











Table 2. Spatially averaged cloud frequency for radiatlvely 
interactive and non-interactive clouds (%). 



(A) 

January 

control 

Land Oceans 

(B) 

Inactive thermal 
minus control 

Land Oceans 

(C) 

Inactive solar 
minus control 

Land Oceans 

Supersaturation clouds 







N. hemisphere 

47.9 

39.9 

1.3 

1.9 

-1.6 

-3.1 

S. hemisphere 

50.2 

40.8 

0.3 

1.5 

0.5 

-2.1 

G1 obal 

48.6 

40.4 

1.1 

1.7 

-0.8 

-2.5 

Convective clouds 







N. hemisphere 

8.9 

12.0 

-1.6 

4.3 

0.6 

-2.7 

S. hemisphere 

23.0 

8.5 

-2.7 

3.1 

-0.4 

-2.0 

G1 obal 

13.5 

10.0 

-1.9 

3.6 

0.4 

-2.3 

All clouds 







N. hemisphere 

50.8 

43.8 

1.4 

3.7 

-1.5 

-4.1 

S. hemisphere 

57.6 

43.8 

1.5 

2.9 

i 

O 

• 

CO 

-3.2 

G1 obal 

53.1 

43.8 

1.4 

3.2 

-1.3 

-3.6 


The model's convective clouds, which represent cumulus clouds in 
nature, generally form less frequently than supersaturation clouds. However, 
the relative changes in their frequency caused by eliminating rari ! ative 
interactions are much larger. The frequency changes over the oceans are cor- 
related in a negative sense with the induced atmospheric temperature change, 
i.e., convective cloud frequencies decrease when the atmosphere is heated by 
enhanced solar radiation, and increase when it is cooled through enhanced 
long wave emission. Changes are on the order of 2-4 percent. 









The opposite is true of convective cloud frequencies over land, 
where frequencies increase by about 0.5 percent when there is stronger heating 
{except in the southern hemisphere) and decrease by 2 percent when there is 
enhanced cooling. Globally, the convective frequencies decrease from January 
control value of 11.2 percent to 9.8 percent with solar transparency and 
increase to 12.9 percent with thermal transparency. Since the convective 
cloud frequencies over the oceans are more sensitive to solar radiation 
changes than those over land, model -generated convective cloudiness apparently 
is more strongly correlated with atmospheric temperature and stability than 
with surface temeprature. On a global basis, the combined frequency of all 
cloud types decreases from the control value of 46.9 percent to 44.1 percent 
without solar interactions and increases to 49.6 percent without thermal 
i nteractions. 

Recently there has been active debate over the relative roles of 
the solar and infrared effects of clouds in the global climate system. In 
particular, would an increase in cloud cover be associated with a surface 
warming due to "greenhouse" type processes, or to a surface cooling due to 
albedo effects? Is there a simultaneous redistribution of clouds that accom- 
panies and possibly compensates for differences in the radiation balance by 
total cloud amount changes? 

On the basis of theoretical models of radiative equilibrium and 
radiative energy balances for specified cloud conditions (e.g., Manabe and 
and Strickler, 1964; Manabe and Wetherald, 1967), it appears that increases 
in the amount of low and middle level clouds cause the surface temperature to 
decrease. But it also appears that increasing the upper level (e.g., cirri- 
forin) cloud amount increases the surface temperature, although this latter 
conclusion is more sensitive to assumptions about cloud optical properties 
and cloud temperature. Thus, the albedo effect tends to dominate when changes 
in low cloud amount occur, while the infrared effect dominates when the 
changes involve upper level clouds. Clearly, any discussion of the radiative 
effect of changing gl obal cloud amounts must consider the possibility of 
changes in the distribution and type of cloud as well, a point stressed by 
Cess and Ramanathan (1978). In fact, Cess (1976) suggests that changes in 
global cloudiness occur in such a fashion that the net radiation balance is 
essentially unchanged. 

The results presented here indicate that some geographical redistri- 
bution of cloudiness does accompany changes in the visible or infrared opacity 
of the atmosphere, and redistribution takes place in the vertical direction 
as well. The increase in globally averaged cloud frequency for the transparent 
infrared case is effected by an increase at virtually all levels in the model, 
except for supersaturation clouds at level 2 and convective clouds at level 5. 
The decrease in total cloudiness for the solar transparent case is caused by 
a decrease of frequencies at most levels, although convective and supersatura- 
tion cloud frequency each increase slightly at one level. There is some 
tendency at most levels for the changes in the cloud frequency caused by 
infrared transparency to compensate partially for the changes due to solar 
transparency, but complete cancellation does not occur. 
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A comparison is also made between the results of the GLAS GCM and 
the radiation balance study of Ellis (1978) based on satellite radiometer 
measurements. For the December-February period, Ellis concluded that the 
cloud effect on the radiation balance, N, which was obtained by subtracting 
cloud-free values from values obtained with actual cloudiness was -12, -42, 
and -27 W nr* for the northern hemisphere, southern hemisphere, and globe, 
respectively. For comparison, the same GCM values obtained here were +9, 
-42, and -16 W nr*, which agree well globally and in the southern hemisphere, 
but not in the northern hemisphere. As with the results from the GCM, Ellis' 
values do not apply to smal 1 changes in cloud amount. However, they are 
valuable because they are measurements of actual atmospheric cloud fields 
that contain the true feedback between cloud formation and cloud radiative 
processes. 
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REVISION OF GLOBAL TOPOGRAPHY ARRAYS 
(D. Randall and J. Abeles) 


In preparation for the testing of a new climate model based on 
the C-grid (Randall et al . , 1980), it was necessary to create a new 
global topography array. We averaged the Scripps 1° x 1° topography data 
set (Gates et al . , 1975) onto a 4° x 5° C-grid. Comparing the result with 
the standard (Tift'S B-grid topography, we found that the new topography 
shows much greater relief over all mountainous regions of the world and 
that the land surface elevations of Greenland, in particular, were much 
higher than in the standard GLAS topography array. 

To study these discrepancies more carefully, we recomputed the 
topography for the 4° x 5° B-grid, so that a direct comparison could be 
made with the standard GLAS topography. As shown in Fig. 1., in the new 
array, the highest elevations of the Himalayas are over 1 km greater, the 
Andes are more than 600 m higher, and there are similar large differences 
in many other parts of the world. Fourier analysis along latitude circles 
shows significant differences in low wave numbers (Fig. 2). At some lati- 
tudes, the topographic forcing at low wave numbers is more than twice as 
strong for the new topography as for the old. 

These findings are of potentially great significance for the GLAS 
climate model i ng effort because, to the extent that the model responds 
linearly to topographic forcing, a doubling of the forcing will produce a 
doubling of the model's response. In fact, one deficiency of past GLAS 
climate simulations has been an underprediction of the amplitudes of the 
long planetary waves. 

Preparations are underway to determine by numerical experiment 
how much impact the revised topography has on climate runs. Preliminary 
tests of the topography in the GLAS 4° x 5° fourth-order forecasting model 
have shown improvements in the forecasts. Revision of the topography used 
in the GLAS 2.5° x 3° forecasting model is also planned. 
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Fig. 1. 


Differences between the land surface height computed from the Scripps topography 
and those traditionally used in the GLAS model. The new topography is higher 


than the old at virtually all land points. 
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TESTING OF A POTENTIAL ENSTROPHY-CONSERVING MOMENTUM 
ADVECTION SCHEME BASED ON THE C-GRID 

(D. A. Randall, J. Abeles, and C-H. Moeng) 


Early attempts to simulate the general circulation of the atmosphere 
(Phillips, 1956) were plagued by nonlinear computational instability, which 
occurred as finite difference truncation errors allowed a false cascade of 
energy to the smallest resolvable scales (Phillips, 1959). The instability 
can be prevented by selectively smoothing the smallest scales, but objections 
can be raised against such smoothing. When the smoothing increases the simu- 
lated rate of dissipation (e.g. Smagorinsky, 1963), it is sometimes interpreted 
as a parameterization of unresolved subgrid scale diffusion, but in practice 
the smoothing techniques are ad hoc and are not derived from turbulence 
theories. On the other hand, if the smoothing is nondissipative, then it 
cannot represent the effects of subgrid scale turbulence. 

Arakawa (1966) showed that the nonlinear instability could be avoided 
through the use of a finite difference momentum advection scheme which conserves 
kinetic energy and, in the limit of two-dimensional nondivergent flow, also 
conserves enstrophy. Arakawa emphasized that his scheme not only prevents the 
nonl inear instabil ity, but also allows simulation of the statistics of the 
general circulation over simulated times far longer than those 'over' which the 
detail s of the circulation can be predicted. In short-to-medium range fore- 
casting (up to about 2 weeks), deterministic prediction of individual synoptic 
features is a credible objective, and one of the most important factors 
determining a model's suitability for such forecasting is the local spatial 
and temporal accuracy of its finite difference scheme. Especially for short 
forecasts, prediction of the statistics of the circulation is not crucial, 
since over such short periods no reasonable scheme will allow the statistics 
to wander far from those of the initial condition. But for climate modeling, 
the required simulation times are so long that deterministic prediction is 
hopeless, and therefore local accuracy is almost irrelevant. Instead, the 
statist ical accuracy of the scheme mostly determines the model's predictive 
ski'lTT 

The present GLAS climate model is based on Arakawa's (1966) scheme, 
and on other aspects <pf ths version of the UCLA model described by Arakawa 
(1972). In spite of the demonstrated success of the present scheme, the 

following three difficulties remain. First, in the limit of two-dimensional 
nondivergent flow, the present scheme conserves enstrophy only in case the 

grid increments are spatially uniform; on a spherical grid, this is true only 

near the equator. Second, as discussed by Arakawa and Lamb (1977), the 

B-grid used in the GLAS climate model is less suitable for simulation of 
the geostrophic adjustment process than the C-grid. Finally, as discussed 
by Arakawa and Lamb (1980), the present scheme does not guarantee conservation 
of potential enstrophy in a divergent flow. As a result, the scheme is not 
wellsuited to simulation of the strongly divergent circulations forced by sharp 
topographic features. 
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Suarez and Arakawa (1979) have presented preliminary results from a 
new version of the UCLA GCM, which incorporates a revised finite difference 
scheme based on the C-grid and designed to conserve potential enstrophy on the 
sphere. The new scheme is thus free from the three deficiencies of the current 
GLAS scheme. But the new UCLA model also incorporates several additional major 
new features, including a vertical coordinate system in which the planetary 
boundary layer top is a coordinate surface. For this reason, it is difficult 
to isolate the impact of the new horizontal differencing scheme on the model 
results. 


We have incorporated the new UCLA horizontal differencing scheme into 
the GLAS climate model, without making any other changes. By comparing the 
results of this new model with those of the current GLAS climate model, we plan 
to identify and evaluate the impact of the new scheme through a series of 
controlled experiments, including both short-range forecasts and climate simu- 
lations. At present (February 1980), all coding changes have been completed, 
and the model has been successfully tested in a 7-day run, starting from 
idealized initial conditions and using idealized topography, without any 
physical parameter! zat ions except dry convective adjustment. 

We have also constructed shallow water models based on the current 
GLAS scheme and on the new scheme. These models incorporate all the essential 
new elements of the respective schemes, but little else. They therefore serve 
as ideal test beds for the schemes, and a large number of experiments have 
already been performed with them. These experiments have shown that the new 
scheme is considerably less noisy than the old, and that even such gross 
statistics as the simulated zonal mean water depth differ significantly 
between the two schemes. Further tests are planned, with varying horizontal 
resolutions, to establish which scheme requires less resolution for a given 
accuracy of approximation to the "exact" (high-resolution) solution. We 
also plan to test a fourth-order version of the new scheme. 
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A NON-ITERATIVE CALCULATION PROCEDURE FOR SURFACE 
TEMPERATURE AND SURFACE FLUXES IN THE GLAS GCM 

(Y. Sud and J. Abeles) 


The present GLAS model has been described by Somerville et al . 
(1974), Stone et al . (1977), and more recently by Hal ern et al. (1979). Its 
surface fluxes of sensible and latent heat were found to exhibit strong 2 -$t 
oscillations at the individual grid points as well as in the zonal and hemi- 
spheric averages. In addition, it was pointed out by Charney et al. (1977) 
that a basic weakness of the GLAS model has been its lower evaporation over 
oceans and higher evaporation over land in a typical monthly simulation. On 
examining the planetary boundary layer (PBL) parameterization, it was found 
that the approximations used in the calculation of surface temperature and the 
use of ad hoc constants in the eddy diffusivity calculation for the mixed 
layer were primarily responsible for these deficiencies. 

By using a quadratic equation for a stable PBL and by a curve-fit 
relation for an unstable PBL, the GLAS model PBL parameterization has been 
changed to calculate the mixed layer temperature gradient. The basic formulae 
used to determine the drag coefficient, its stability dependence and the effect 
of moisture on the temperature gradient remain unchanged. This calculation in 
the PBL removed the 2-6t oscillations in the surface temperature and surface 
fluxes. Also, the geographical distributions of the surface fluxes were 
improved. 


Another advantage of this calculation was to avoid iteration in the 
surface flux calculation thereby producing some speed improvements. Compari- 
sons of the evaporation over land and ocean are shown for the old and new 
calculations. The ba|ic quantities and useful variables are defined below: 

Tg - ground temperature, °K 

qg - saturation mixing ratio at the ground 

T 9 , 9 = physical and potential temperature of layer 9 in the model 
qg = mixing ratio at level 9 of the atmosphere 
q s = surface mixing ratio 
U $j V s = surface U,V winds, in ms"^ 

Ug Vg = level 9 U,V winds, in ms“^ 

W s = (U s 2 + V $ 2 )V2 

K - Eddy diffusivity of the mixed layer cal nHk " 1 
L = latent heat of evaporation calg'l 


169 


a 








'viiiiW***' a. 1 - - 


Cp * specific heat at constant pressure 
R * gas constant 


k, * R/Cp 

e * ratio of molecular weights of water and air 


Fig. 1 shows a schematic sketch of a potential temperature profile in a stable 
and an unstable PBL, The temperature difference across the mixed layer, , 
is modified to reflect the moist adiabatic lapse rate as follows; 


Y * Yd + ^s SL aZ B 1 
c c p 


* 



1 + 5418 


Define 


and 


A T^ B F p » AZg , 

s 

a Td * r d a Zb i 


Y c s r s .2. Zg 
C P 


L 

t - 1 * cpir 



L 

1 + 5418 ^ 



t 


The above equations yield 


Tc = Td + Yc * 

Assuming that Dr is prescribed on a grid cell, the surface heat 
flux, F*, may be obtained as follows: 

Fs= pCpDR(0g- 0 S ) 


where 

d r = c d w s 3 /( w s 2 - 7x(Ae- 8)) 
Dr = Co(W s + /A0^ T ) 


if T g <T s 
if T g >T s . 
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The mixed-layer heat flux, F m i , is given by 


F m l a - Cp K/Z B (AT -A T c ) , 


First, for the stable case 


K * 60. /(I + 40.*Ri) , 


where 

r. = -6A gZ B /[0 g * (U 9 -U s ) 2 + (Vg - V s ) 2 ]. 
Second, for the unstable case 


K = 60. + 100. rl - exp(-l .2 60. )[=] 

5z 

k By equating F s and F [n , we find that 

i 


D R (8g - 0 s ) = -K/Z B ( AT - AT d - y c ) « K/Z B (d 0 + Yc )- 
( a ) S olution of 69 (Stable Case) 

The following equation for the stable case is obtained by using the 
above equations: 


i 

! 


r 


/ 


(7 x 60/Z B C D W s 3 - 40gZ B / g D w 2 )60 2 + 

^ Y - 

A 

(1 + 60/Z b C d W s - 7x60 /Z B C D W S 3 + 40gz B C D / gD w 2 )^0 =A0 , 

V y * ' ' 

B 


from which 59 may be obtained. However, if A approaches zero in the above 
quadratic equation, the limiting solution is 


6 0 = -A0/B(1 - A.A3/B 2 ). 


The boundary layer height Z B 
stable PBL is generally shallow. Its 
limiting value of K = 2.0 is reached 


is of the order of 500 m. However, a 
height is of the order of 100 m. A 
when the limiting critical Richardson 
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number 3.05 is attained (Deardorff, 1967). Also, the constant (40.) is 
adjusted to (8.) to reflect use of 500 m for the PBL height. The minimum 
value of 60 is obtained by 

se =Ae/n. + (wWsio - 7- A ?« s 2 n • 


(b) Solution of (Unstable Case) 

In an unstable case, a transcendental is involved. Therefore, an 
exact solution is not possible. Hence, an attempt is made to obtain 59 as a 
function .of A9 , w s , and Cq, for a range of values. Figs. 2a, 2b, and 2c show 
6 8 as a function of A9 for values of Cq between .001 and .005 and surface wind 
magnitudes of 2-12 m s - ^ . From these solutions, graphs are obtained for 6 8 / A 0 
versus C^, for various values of the winds (Fig. 2d) and 69/ A0 versus surface 
wind for various values of Cq. A relation between 60 andA9 , Cp and W s , was 
curve fitted to give 


60 = (0.1382 . W s . C D 1/2 + 13.67 C D ) 


(c) Wind Field Modification by Surface Drag 

In the old code, first the wind fields were interpolated to primary 
points, then momentum deficit was calculated for the interpolated wind fields. 
Having done that, the momentum deficit was reinterpolated to secondary points. 
This back-and-forth interpolation does not allow a direct coupling between 
wind velocity and momentum deficit. Replacement of this procedure by the one 
used in the new boundary layer calculation helps to partially eliminate a 2- x 
oscillation. In the new procedure, first a factor, , is calculated at each 
primary point, as before. It is then interpolated linearly to secondary points. 
The momentum deficit is now calculated by multiplying this factor by winds at 
level 9. Also, wind velocity updates are saved in an array. These are made 
simultaneously at the end. 

Analysis of Results 

(1) First, a time step invoking one call each of radiation and 
physics and three calls of hydrodynamics was completed to compare the results. 
The only striking difference was larger evaporative tendencies over oceans and 
some large differences in surface temperature. Fig. 3 shows the digitized 
maps of surface temperature differences. A large difference of 10-20°C may be 
seen in regions marked (I). These differences occur because in this region 
the air aloft the surface is very cold. The old model makes T s close to the 
ground temperature, whereas the new model makes it close to level 9 potential 
temperature. This is so because in the new calculation the eddy exchange 
coefficient has increased by an order of magnitude for the unstable atmosphere. 

In a 1-day simulation, the 2-§t oscillation in the surface temperature 
and surface fluxes was eliminated over land and ocean. See Figs. 4a and 4b 
for changes in the fluxes of evaporation and sensible heat over land and 
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ocean in the two hemispheres separately. Tho previous PBl model runs exhibit 
a very noisy field compared to the new PBL model. Consistent changes in daily 
averages of sensible heat and evaporative fluxes on individual grid points are 
now simulated. 

The new PBL parameterization has been tested with a unique solution 
of surface temperature instead of an iterative and time-averaged solution. 
Also, the original formulation of the eddy diffusion coefficient of Deardorff 
(1967) has been incorporated. The results from the new parameterization show 
desirable effects of increased evaporation over ocean, reduced evaporation 
over land along with elimination of 2-St oscillation in the surface fluxes of 
the old model. Also, a short-range forecast revealed small but beneficial 
effects on surface temperatures, sea level pressure and geopotential heights 
at 500 mb. 
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Fig. 3. Differences between surface temperature calculated at the initial 

time from the old and new formulation of PBL transport calculations. 
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SPACE-TIME SPECTRA OF THE GLAS MODEL 500 MB GEOPOTENTIAL HEIGHT 
(D. M. Straus and J. Shukla) 



I 


Space-time spectral analysis is a technique which expresses the 
space-time variance or covariance of meteorological fields in terms of con- 
tributions from individual zonal wavenumbers and frequency bands. Because 
motions of distinct space and time scales are often associated with different 
dynamical origins, space-time spectral analysis provides a useful diagnostic 
tool with which to compare model output with observations (Hayashi and Golder, 
1977; Pratt, 1979). 

This paper represents part of a continuing effort to document the 
winter space-time spectra of the GLAS general circulation model (Halem et al ■ , 
1978) and to compare them with the spectra of 14 winter seasons obtained from 
the NMC analyses for the years 1963-64 through 1976-77. In this paper, only 
the 500 mb geopotential height is considered. 

The spectra can best be introduced by consideration of the quantity 
G(A,t), the 50Q mb geopotential height as a function of longitude (x), and time 
(t) for a fixed latitude. It must be understood that the seasonal mean and 
some estimate of the seasonal cycle has already been removed. If this pro- 
k cedure is not followed, the well known phenomenon of "leakage" will cause a 

; gross distortion of the resulting seasonal spectra. (For a brief discussion 

: of the phenomenon of leakage, see Bingham et al . , 1967.) The space-time vari- 

ance of G is defined as 

1 

»*if /-‘M a . in 

o o 

i where T is the period under consideration (here T * 90 days). V is further 

written as 

00 

V - £ I P(k.f) • (2) 

k=0 f>0 

Here k is the (discrete) zonal wavenumber and is constrained to be nonnegative, 
while the (discrete) frequencies f are given by 

f “ £ > 1 < m < 90 , (3) 

T 

for twice daily data. The variance P(k,f) can be further written (for k^O) as 

P(k,f) = PW(k,f) + PE(k,f ) , (4) 

where PW is the valance due to westward propagating motions and PE is due to 
eastward propagating motions. The meaning of PW and PE can best be appreciated 
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by noting that any function of longitude and time, and in particular G, can be 
written as 


« ]T Z (A w Oc,f)cos(kX+2Ttft+(j) w ) 

k«o f o 


( 5 ) 


+A E (k,f)cos(kX- 27Tft+<}> E ) j , 


where Ay and A E are westward and eastward propagating waves, respectively. 
Ay and A E are then related to PW and PE of Eq. (4) by 

PW(k,f) ~ **Ay(k,f) 

2 ( 6 ) 
PE(k,f) * >iA E (k,f) . 


The tendency for either westward or eastward motions to dominate the total 
variance is given by the "net propagation tendency" 


PR(k.f) 



(7) 


It is of import a rice to clearly delineate what kind of information is 
obtained from PW and PE. They do not, by themselves, indicate whether the 
fluctuations are due to physical standing waves, propagating waves, or neither. 
They only reflect the projection of the total (geopotential) field onto a 
complete set of basis functions which correspond to traveling waves. Thus, 
for instance, a value of PR close to -1.0 indicates a strong tendency f c " 
motions of a particular space-time scale to propagate eastward, but a value 
near 0.0 only means that no clearly dominant direction of propagation is 
apparent. 

The model data were obtained from a winter simulation of the GLAS 
model which was initialized with data from 1 January 1975 and consist of twice- 
daily geopotential heights for the first 90 days of the simulation. The NMC 
analyses were obtained from NCAR (see Jenne, 1975) and consist of twice-daily 
fields obtained for 14 winter seasons (1963-64 through 1976-77). Each winter 
season was considered to start on 1 December and to be 90 days in length 
In order to approximate th> seasonal cycle, a parabolic was fit to each time 
series (whether from the model or the NMC analyses). This parabola was 
removed, and the remaining time series was tapered by applying a half cosine 
bell to the first and last 10 percent of the series. This procedure, suggested 
by Bingham et al . (1967), insures that those components of the annual cycle 
which were not accounted for by the (admittedly crude) parabolic approximation 
would cause minimal distortion of the calculated spectra. 

The space-time variance P(k,f), obtained from the NMC analyses, is 
presented in Fig. 1. The function shown results from averaging the functions 
P calculated for each of the 14 seasons separately. Two important features of 
P must be noted. One is the concentration of variance in the lowest frequency, 
lowest wavenumber portion of the spectrum. The other is the presence of a 
"ridge" of high variance extending from low frequencies and wavenumbers to 
higher frequencies and wavenumbers. In Fig. 1, the dotted lines delineate 





four space-time regions which divide the wavenumber-frequency plane into 
regimes of (i) low frequency planetary waves; (ii) low frequency medium-scale 
waves; (iii) higher frequency planetary waves; ar.J (iv) higher frequency 
medium-scale waves or baroclinic-scale waves. The precise definitions of 
these regions are given in Table 1. Region (i) has the bulk of the variance 
of the planetary waves, and region (iii) consists of space-time scales roughly 
corresponding to the most unstable baroclinic motions. The remainder of this 
paper will concentrate on these two scales. 


Table 1. Definitions of spectral regions. 


Region 

Range of 
frequencies 
(cpd) 

Range of 
period 
(days) 

Range of 
wavenumbers 

(i) 

.0111 to .1333 

90 to 7.5 

1 to 4 

(ii) 

.0111 to .1333 

90 to 7.5 

5 to 10 

(iii) 

.1444 to .3556 

6.9 to 2.8 

1 to 4 

(iv) 

.1444 to .3556 

6.9 to 2.8 

5 to 10 


The total space-time variance in region (i) is presented for. the 
northern hemisphere in Fig. 2, both for the GLAS model and the 14-winter NMC 
average. The number of degrees of freedom for all of the calculations are 
summarized in Table 2. It is clear that the model underestimates the average 
low frequency planetary wave variance in midlatitudes by a considerable amount. 
In addition, the observed maximum between 60N and 70N is absent in the model. 
In Fig. 3, the model's low frequency planetary wave variance is compared to 
that of two individual winters. These winters were chosen to define a range 
of interannual variation of the variance in region (i); the winter of 1967-68 
was a relatively active one in terms of low frequency planetary waves, while 
the winter of 1969-70 was a relatively quiet season. The model underestimate 
of the region (i) variance is marked, even when comparing to the latter season. 

The baroclinic-scale variance (region [iii]) of the model is compared 
to the averaged observations in Fig. 4. The agreement is excellent, both as 
to the overall magnitude of the variance and its dependence on latitude. In 
Fig. 5, the comparison is again made between the model variance and that of 
two individual winter seasons chosen to define the range of interannual varia- 
bility. The model variance is much more similar to that of the barocl inical ly 
active winter than it is to the quieter one. The propagation tendency PR, 
averaged over the baroclinic space and time scales (region [iii]) is presented 
in Fig. 6 for the model and the same two winter seasons. The model agrees 
well with both seasons in that the propagation tendency is close to -1.0, 
indicating a strong tendency to propagate eastward. 
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Fig. 2. Variance in regime (i) (low frequency planetary waves) for 
the wdel (dots) and for the average of 14 winter seasons 
(sol id 1 1 ne) . 
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Fig. 3. Variance in regime (i) (low frequency planetary waves) for 
the model (dots), for observations from the winter of 1969-70 
(dashed line), and for observations from the winter of 1967-68 
(solid line). 
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Fig. 4. Variance in regime (iv) (barocl i nic-scale waves) for the model 
(dots) and for the average of 14 winter seasons (solid line). 
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1965-66 

1974-75 


LATITUDE 


Variance in regime ( 1 v ) (barocl inic-scale waves) for the model 
(dots), for observations from the winter o* 1974-75 (dashed 
line), and for observations from the winter of 1965-66 (solid 
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Fig. 6. Average propagation tendency In regm* ( 1 v ) (barocl Inlc-scale 
waves). As In Fig. S. 
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Two salient points emerge from the results that have been presented. 
One is that the GLAS model simulates baroclinie-scale disturbances in a very 
realistic manner, at least in a statistical sense. The second point is that 
the model seriously underestimates the amount of variance in the low frequency 
planetary waves. Comparisons of the GFDL and NCAR models with observations 
show that these models share the latter deficiency (underestimate of region 
[i] variance)-, however, they have not had the GLAS model's success in repro- 
ducing baroclinic-scale waves (Pratt, 1979). 


Table 2. Degrees of freedom. For a single season, the number of degrees of 
freedom is 2 per wavenumber per harmonic frequency. The effect of 
the cosine bell taper is to reduce the number by a factor of about 
0.87 (Julian, 1971). 


Region 

Model 

Observations 
(single season) 

Observations 
(14-year average) 

(i) 

84 

84 

1169 

(ii) 

125 

125 

1754 

(TIT) 

139 

139 

1949 

(iv) 

209 

209 

2923 


The spectra that have been presented describe the transient motions 
of the atmosphere and the model, exclusive of the annual cycle. It is also of 
some interest to study the mean fields. These were computed by averaging the 
(unadjusted) model or observed fields over the appropriate 90-day season. The 
total spatial variance of the mean field can be written as a sum of contribu- 
tions from each zonal wavenumber. The sum of the variances in wavenumbers 1 
through 18 (stationary waves) is presented in Figs. 7 and 8. In Fig. 7, the 
model is compared to the average of 14 seasons, while in Fig. 8 the model is 
compared to two winter seasons, one in which the stationary waves are strong 
(1976-77) and one in which they are weak (1970-71). The results show that the 
model's stationary waves are far too weak in midlatitudes. This is in contrast 
to the GFDL and NCAR models, whose simulated stationary waves are in better 
agreement with the observations (Pratt, 1979). 

In conclusion, the space-time spectral analysis of the 500 mb height 
field has indicated some important successes and failures of the GLAS model 
in simulating the atmosphere. The model does simulate the baroclinic-scale 
waves well, both in terms of the total variance and the dominant direction of 
propagation. On the other hand, the model's stationary- and planetary-wave 
variances are much too small compared to the atmosphere. 
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Fig. 7. Stationary-wave variance for the model (dots) and for the 
average of 14 winter seasons (solid line). 
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Stationary-wave variance for the model (dots), for observations 
fro^ the winter of 1970-71 (dashed line), and for ooservations 
from the winter of 1976-77 (solid line). 
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OROGRAPHICALLY INDUCED MULTIPLE EQUILIBRIA IN A 
HIGHLY TRUNCATED BAROCLINIC SPECTRAL MODEL 

(J. G. Charney and D. M. Straus) 


In a recent study, Charney and DeVore (1979) showed that the forced 
flow of a barotropic fluid over wavy topography In a periodic beta-plane channel 
may possess a multiplicity of equilibrium states of which more than one may be 
stable. In a certain range of forcing, there are two stable and one unstable 
equilibrium states. One of the stable states is a "high-index" circulation 
having a strong zonal flow and a relatively weak wave perturbation; the other 
is a low-index flow having a weak zonal flow and a relatively strong perturba- 
tion, i.e., the unstable equilibrium state has an intermediate zonal flow. It 
was suggested in Charney and DeVore that blocking, the presence of persistent 
high-pressure ridges in the Pacific, Atlantic or Central Asiatic regions is 
an example of an equilibrium flow of the low-index type, and that transition 
from it to the high- index state occurs via the instability of the intermediate 
state. 


The discovery of multiple equilibria in a barotropic system led the 
present authors to investigate a corresponding baroclinic system. Are there 
multiple stationary equilibria, and, if so, what is their nature? Holopainen 
(1970) showed that the quasi-stationary planetary waves in the atmosphere are 
maintained by conversion of mean flow potential energy by baroclinic processes, 
not by conversion of mean flow kinetic energy by the action of mountain pressure 
torques. What then is the role of orography? 

To investigate these questions, the authors studied a highly trun- 
cated, two-layer, quasi-geostrophic spectral model of the type introduced by 
Lorenz (1963). The model domain is a periodic channel on a mid-latitude beta 
plane. The spectral equations studied are precisely those of Lorenz, with the 
addition of topography and a variable Coriolis parameter, and the omission of 
an equation governing the time variation of the static stability. The param- 
eters (such as the channel width, frictional dissipation, static stability, 
etc.) were taken to be representative of the real atmosphere. The modes re- 
tained correspond to a zonal flow plus only one zonal wave number, and to two 
modes in the north-south (y) direction. (See Charney and Straus, 1980, for 
further details. ) The results discussed in this report are those for which 
the single zonal wave number was taken to be zonal wave number 3 in the real 
atmosphere at latitude 45N. The orography was taken to consist of a single 
(wave number 3) sine wave in the x-direction, and to have a projection onto 
only the gravest (first) y-mode. 

The Hadley solution is an equilibrium state obtained by setting all 
time derivatives and wave amplitudes to zero. The zonal flow in the lower 
layer also vanishes, so that the Hadley solution remains a true equilibrium 
in the presence of orography. 

An analysis of the stability of the Hadley solution to (wave number 3) 
perturbations of the first y-mode reveals the presence of two types of insta- 
bility. They are distinguished by the character of the eigenvalue of the 
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stability matrix whose real part becomes positive, i.e., the eigenvalue cor- 
responding to the instability, In one case, this eigenvalue is purely real, 
corresponding to a perturbation that grows in place. The instability is termed 
"orographic." In the second case, the eigenvalues occur in complex conjugate 
pairs, indicating a traveling disturbance. This Instability is a modified 
form of baroclinlc instability. For the set of (realistic) parameters used by 
Charney and Straus (1980), the Hadley solution is stable for thermal driving 
O^less than .131. The parameter 0 A * is a dimensionless measure of the radia- 
tive equilibrium temperature difference between the southern and northern 
walls. It is related to the dimensional radiative equilibrium temperature 
difference (AT*) and to the actual temperature difference (AT) that occurs in 
the Hadley solution, by 


AT* - 530 0 A * (°K), 

AT * 3780a* (° k )‘ 

For 0/\*>,131 the Hadley solution becomes orographical ly unstable 
(see Table 1), and for 8 A * above ,165 it also becomes baroclinically unstable. 
Thus, in the range of forcing . 1 65< 0 A *<,2OO, the two types of instability 
coexist. 


For values of o A * below .131, the lower limit of topographic insta- 
bility, no equilibria other than the Hadley circulation exist. At values of 
the driving which are slightly supercritical (0 A * « .135), two new equilibria 
are found, These are referred to as "wavy equilibria" because they are steady 
state solutions with non-zero wave amplitudes. (For details regarding their 
computation, see Charney and Straus, 1980.) 

Each solution is identified by a branch number. (See Tables 1 and 2.) 
The two new equilibria are themselves stable to further first y-mode perturba- 
tions. Characterizing the three steady states by the zonal flow in the lower 
layer, the (unstable) Hadley solution (zero lower layer zonal flow) lies inter- 
mediate between the (stable) branch 1 solution (weak lower layer westerlies) 
and the (stable) branch 2 solution (weak easterlies). Since both new solutions 
are stable, they divide the phase space corresponding to the first y-mode into 
two subdomains. 

Numerical integration of the equations of motion (with only the first 
y-mode present), starting near either equilibrium, shows convergence to that 
equilibrium. Both "wavy" equilibria are, however, unstable with respect to 
second y-mode perturbations. Numerical integration with both y-modes present 
led to a variety of solutions, depending upon the initial conditions. It was 
possible to obtain a periodic solution (with a period of about 16 days) which 
consisted of a relatively weak stationary wave, upon which was superimposed 
transient waves. The first y-mode transient wave propagated westward whereas 
the second y-mode wave propagated eastward. 

As the value of e A * is increased, the two wavy equilibria become 
modified in a continuous manner. Branch 1 is characterized by lower layer 
zonal westerlies and branch 2 by lower layer zonal easterlies, characterizations 
which remain valid throughout the, range of 8 A * studied. 


V-U 
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Table 1. Wavy equilibria for wave number 3, branches 1 and 2. 






Branch 1 





s? 5 

* 

Sh 


% 


0 A 

0 K 

1 

01 

S 1 

S2 

.135 

Or 

.0982 

-.0090 

.0174 

.0919 

-.0055 

.0169 


P 

.140 

Or 

.1000 

-.0167 

.0231 

.0908 

-.0111 

.0223 


P 

.145 

Or 

.1014 

-.0237 

.0265 

.0899 

-.0163 

.0255 

-- 

P 

.150 

Or 

.1028 

-.0302 

.0287 

.0891 

-.0210 

.0274 

-- 

P 

.155 

Or 

.1041 

-.0362 

.0300 

.0885 

-.0253 

.0286 

P 

P 

.160 

Or 

.1055 

-.0421 

.0309 

.0879 

-.0293 

.0294 

P 

P 

.165 

Or, Be 

.1067 

-.0475 

.0316 

.0874 

-.0328 

.0299 

P 

P 

.170 

Or, Be 

.1078 

-.0527 

.0320 

.0870 

-.0361 

.0302 

P 

P 

.175 

Or, Be 

.1090 

-.0577 

.0323 

.0866 

-.0392 

.0303 

P 

P 

.180 

Or ,Bc 

.1100 

-.0622 

.0326 

.0863 

-.0419 

.0304 

-- 

P 

.185 

Or, Be 

.1111 

-.0667 

.0328 

.0860 

-.0445 

.0305 

-- 

P 

.190 

Or, Be 

.1120 

-.0709 

.0329 

.0857 

-.0468 

.0306 

-- 

P 

.195 

Or, Be 

.1129 

-.0751 

.0331 

.0855 

-.0491 

.0306 


P 

.200 

Or ,Bc 

.1133 

-.0/70 

.0334 

.0855 

-.0500 

.0309 

-- 

P 





Branch 2 







1 





1 



0 A * 

Sh 

*A 

'K 


°A 

e K 

0 L 

Si 

s 2 

.135 

Or 

.0930 

-.0001 

-.0057 

.0960 

-.0010 

-.0055 


P 

.140 

Or 

.0927 

-.0024 

-.0125 

.09/7 

-.0046 

-.0120 

P 

P 

.145 

Or 

.0918 

-.0055 

-.0164 

.0993 

-.0086 

-.0157 

P 

P 

.150 

Or 

.0910 

-.0088 

-.0186 

.1008 

-.0126 

-.0177 

P 

P 

.155 

Or 

.0904 

-.0120 

-.0197 

.1026 

-.0165 

-.0186 

P 

P 

.160 

Or 

.0899 

-.0150 

-.0198 

.1047 

-.0201 

-.0185 

P 

P 

.165 

Or, Be 

.0895 

-.0177 

-.0192 

.1071 

-.0232 

-.0176 

P 

P 

.170 

Or, Be 

.0895 

-.0198 

-.0179 

.1100 

-.0257 

-.0161 

P 

P 

.175 

Or, Be 

.0898 

-.0214 

-.0164 

.1135 

-.0274 

-.0142 

P 

P 

.180 

Or, Be 

.0907 

-.0224 

-.0146 

.1174 

-.0283 

-.0122 

P 

P 

.185 

Or, Be 

.0924 

-.0229 

-.0128 

.1217 

-.0284 

-.0102 

P 

P 

.190 

Or, Be 

.0949 

-.0230 

-.0110 

.1264 

-.0277 

-.0082 

P 

P 

.195 

Or, Be 

.0990 

-.0225 

-.0092 

.1315 

-.0260 

-.0063 

P 

P 

.200 

Or, Be 

.1052 

-.0212 

-.0071 

.1370 

-.0229 

-.0043 

s,s 

P 


All quantities are in dimensionless form. The column headed describes the 
stability of the Hadley circulation to perturbations of the first y-mode. The 
columns S] and S2 describe the stability of the wavy equilibria to perturba- 
tions of the 1st and 2nd y-modes, respectively. No entry (--) denotes a stable 
equilibrium, Or denotes orographic instability of the Hadley circulation and 
Be baroclinic instability of the Hadley circulation. P denotes a propagating 
instability of the wavy equilibrium (complex conjugate roots of the stability 
matrix), S denotes a stationary instability (purely real roots). When two types 
of instability are indicated, they exist simultaneously. The meaning of the 
variables follows: the lower layer zonal flow is T A - 0 A ; the upper layer zonal 
flow is T A + 9/\. The wave in the lower layer is proportional to (fK-0K) cos kx + 
(TL - \)sin kx, where k is the wave number. The wave in the upper layer is 
proportional to (T«; + 0 K )cos kx + (\ + ®l) sin kx. 
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Table 2. Wavy equilibria for wave number 3, branches 1 and 4. 


V 

S H 



Branch 3 



Si 

f~ 

1 # 
1 


\ 

\ 

0 A 


«L ' 

.180 

Or, Be 

.1602 

.0633 

.0024 

.1134 

.0327 

-.0017 

s 

p 

.185 

Or, Be 

.1608 

.0802 

.0019 

.1087 

.0398 

-.0027 

s 

p 

.190 

Or, Be 

.1606 

.0902 

.0016 

.1064 

.0441 

-.0033 

s 

p 

.195 

Or, Be 

.1609 

.0995 

.0009 

.1045 

.0477 

-.0040 

s 

p 

.200 

Or, Be 

.1603 

.1051 

.0008 

.1035 

.0503 

-.0043 

s 

p 





Branch 4 





0 A * 

$H 

r 

*A 


\ 

°A 

®K 

1 

e L 

Si 

S2 

.180 

Or, Be 

.1573 

.0383 

.0256 

.1224 

.0217 

-.0006 

p 

p 

.185 

Or, Be 

.1524 

.0212 

.0020 

.1300 

.0132 

-.0000 

p 

p 

.190 

Or ,Bc 

.1471 

.0103 

.0012 

.1352 

.0069 

.0001 

p 

p 

.195 

Or, Be 

.1411 

.0014 

.0002 

.1393 

.0010 

.0000 

p 

p 

.200 

Or, Be 

.1324 

-.0073 

-.0013 

.1425 

-.0059 

-.0004 

s,p 

D 


All quantities are in dimensionless form. The column headed describes the 
stability of the Hadley circulation to perturbations of the first y-mode. The 
columns S] and S2 describe the stability of the wavy equilibria to perturbations 
of the 1st and 2nd y-modes, respectively. No entry (--) denotes a stable 

equilibrium, Or denotes orographic instability of the Hadley circulation and Be 
baroclinic instability of the Hadley circulation. P denotes a propagating 
instability of the wavy equilibrium (complex conjugate roots of the stability 
matrix), S denotes a stationary instability (purely real roots). When two types 
of instability are indicated, they exist simultaneously. The meaning of the 
variables follows: the lower layer zonal flow is \ - 0 A , the upper layer 
zonal flow is ¥/\ + 6/\. The wave in the lower layer is proportional to 
(\ - 0 k) cos ^x + Of l " 0|J sin k x » where k is the wave number. The wave in 
the upper layer is proportional to (¥« + 9«)cos kx + (4^ + 0|J sin kx. 

At 9/\* = *140, the branch 2 equilibrium becomes unstable with respect 
to first y-mode perturbations. Since the roots of the stability matrix occur 
in complex conjugate pairs, the disturbance has a traveling nature, and in fact 
can be considered as a superposition of a westward and an eastward traveling 
wave. Thus, the presence of topography allows for instability of the wavy 
equi 1 i bri urn even though the Hadley solution itself is barocl inically stable . 
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For still higher values of the driving ( 9 a* = *155), both wavy 
equilibria become unstable to first y-mode perturbations. However, the branch 1 
solution, which by now has acquired a larger upper layer wave, is only weakly 
unstable (with e-folding times of about 200 days). Numerical integration of 
the first y-mode equations gives rise to a periodic solution (with a period of 
about 17 days) in which the transient wave again propagates westward. 

At = .165, the Hadley circulation becomes both baroclinically and 
orographically unstable and, at « .180, two new wavy equilibria are found, 
giving four branches in all. (See Table 2.) Evidence is presented in Charney 
and Straus (1979) which indicates that the existence of both baroclinic and 
orographic instability of the Hadley circulation is closely connected to the 
appearance of new wavy equilibria. These new equilibria can be characterized 
as high index, having very large upper layer zonal flows (corresponding to 
about 60 m/sec) and relatively much weaker wave fields. Both equilibria have 
zonal westerlies in the lower layer. In contrast, the branch 1 equilibrium 
evolves into a blocking, or low-index solution, with an upper layer wave 
amplitude comparable to the upper layer zonal flow. 

The high-index solutions (branches 3 and 4) as well as the branch 2- 
solution are unstable with respect to perturbations of the first y-mode. The 
branch 1 solution, which is weakly unstable for .155< ^*_<.175, becomes stable 
again for 0 a*>. 18Q. Thus, for the large value of the thermal driving cor- 
responding to"V = *180, the only realizable steady state solution of the 
first y-mode system is the low-index, or blocking solution. Numerical integra- 
tion of the nonlinear equations, with initial conditions close to either the 
branch 1 or branch 2-equilibrium, converges to the stable low-index equilibrium. 
When the initial conditions are close to the branch 3 equilibrium, a periodic 
solution of 18.1 days is realized, and when they are close to the branch 4 
steady state, a similar but not identical periodic solution (of period 16.6 
days) is obtained. The time-averaged zonal flow in both cases is similar to 
that of the stable equilibrium, but the time-averaged wave amplitudes are 
weaker. The time-dependent flow consists of slow, westward traveling waves of 
variable amplitude accompanied by large oscillations in the zonal flow. All 
four equilibria are unstable to second y-mode perturbations, and both periodic 
and aperiodic solutions are possible when both y-modes are retained, depending 
upon the initial conditions. 

In summary, orographic instability of the Hadley circulation to 
disturbances of wave number 3 gives rise to a wavy equilibrium which for 
large driving resembles blocking. However, baroclinic instability of the 
Hadley solution is necessary to realize high-index equilibria. 

To understand better the role of the mountains in the formation of 
the wavy equilibria, the energetics of their maintenance has been examined. 
The important result is that the orographic form drag exchange between zonal 
and eddy kinetic energies (the only such exchange for the first y-mode) is in 
all cases small. For lower values of the driving, this exchange is entirely 
negligible, while for high driving the small orographic exchange always draws 
energy from the wave. Thus, the orographic interactions play an indirect role 
in the formation of wavy equilibria; they are a catalyst; they permit the 
equilibria to exist, but they do not directly support the stationary waves. 
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Although it is not possible to arrive at definite conclusions con- 
cerning the real atmosphere or even a continuous model from the highly truncated 
spectral, two-layer model examined, a number of heuristic conclusions are 
clearly suggested. One conclusion is that the phenomenon of blocking is a 
stable or quasi-stable equilibrium arising from orographic instability with 
strong thermal driving; and that high-index states occur also for strong driving, 
when the Hadley circulation becomes baroclinically as well as orographically 
unstable. These latter will always be unstable to large-scale perturbations, 
and will give rise to a wave-zonal flow mean state on which is superimposed 
propagating waves, drawing on the potential and kinetic energies of the mean 
flow and interacting with it in a periodic or at least fluctuating manner. 
The results further suggest that the westward propagating first y-mode waves 
arising from the instability of the w, vr- : » <e equilibria are prototypes of low 
frequency, westward propagating "Rossby" wves observed in nature (Eliasen and 
Machenhauer, 1965, 1969, Pratt and Wallace, 1976; Blackmon, 1976; Madden, 1978). 
These waves, which have periods ranging from 5 to 20 days, have been previously 
regarded as "free" waves (Eliasen and Machenhauer, 1965, 1969; Madden, 1978), 
this study suggests that they arise from instabilities of the quasi-stationary 
topographically forced equilibria. 
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A STUDY OF THE CLIMATIC NOISE OF THE 6LAS GCM 
(D. M. Straus and M. Halem) 


The problem of distinguishing the potentially predictable component 
of observed monthly and seasonal atmospheric means from the unpredictable 
statistical fluctuations has recently received much attention (Leith, 1973, 
1975; Madden, 1976; and Madden and Shea, 1978). The potentially predictable 
component (Leith, 1975) is defined as that element of long-term atmospheric 
states which is forced by influences external to the atmosphere, such,as solar 
insolation, sea surface temperatures, ice cover, etc. The unpredictable sta- 
tistical fluctuations (or "climatic noise") are due to the variability intrinsic 
to the atmosphere and would exist even in the presence of unchanging external 
conditions. That is, the statistical properties of climatic noise are defined 
by averaging over an ensemble of realizations, each in the presence of identical 
external conditions. In these definitions, the external influences may be 
fixed, or they may be prescribed to vary smoothly in time (as does the annual 
cycle) . 


However, this distinction is difficult to apply in practice, since 
the operational definition of climatic means is based on finite time averages 
of a single realization, namely, the atmosphere itself. The determination of 
the character of climatic noise from long records of real data requires that 
a priori assumptions be made. One such assumption is that climatic noise can 
be represented as a first-order autoregressive (Markov) process over a wide 
range of frequencies (Leith, 1973). 

In this paper, stochastic forecasts were made with the GLAS GCM to 
test the adequacy of low-order autoregressive processes in describing the 
model climatic noise. Model integrations were carried out with prescribed 
climatological boundary conditions for differing initial states, thus providing 
an ensemble of simulated data which satisfies the definition of climatic noise. 
The fields examined were the sea level pressure and ground temperature at 
54 surface stations in the continental U.S. The hypotheses tested were that 
these data are adequately described by white noise (H 0 ), by a first-order 
autoregressive process (Hi), or by a second-order process (H 2 ). 

The findings show that H 0 and H] must be rejected at virtually all 
stations, but that hypothesis H 2 is consistent with the model "climatic 
noise" in some geographical regions. There are regions for which none of the 
processes considered describe the data. 

The eleven integrations of the GLAS model used in this study were 
performed with the constant external forcings set equal to their January 
climatological values. The only exception is the incoming solar radiation, 
which was also the same in all integrations but changed in time, following the 
declination of the sun. The initial atmospheric conditions for the first of 
the integrations (run 1) were taken from an objective analysis of observed 
data valid for 1 January 1973; four additional sets of initial conditions 
(runs 2-5) were obtained from this analysis by introducing random errors with 
a Gaussian distribution into all fields Similarly, the initial conditions 
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for three other integrations (runs 6-8) were those observed for 1 January 1974, 
plus two sets of similarly perturbed fields. The initial conditions for the 
remaining three integrations (runs 9-11) were those observed for 1 January 1975, 
again with two sets of similarly perturbed fields. Numerical integrations from 
1 January to 31 January were carried out for each of the eleven initial states. 
The hourly sea level pressures and ground temperatures were interpolated to 
the positions of 54 selected weather stations over the U.S. (see Fig. 1) and 
then averaged each day. Thus, for each station and for both variables, the 
data consisted of eleven 31-day time series. In processing these data, it is 
necessary to remove the mean from each of the series separately. Further, it 
is necessary to remove the seasonal trend. An approximation to this trend was 
obtained by averaging the temperature or sea level pressure for each day across 
the eleven ensemble members. The resulting (average) time series was subtracted 
from each of the eleven original series. 

The first hypothesis H 0 is that the data, consisting of N s « 11 time 
series of N = 31 days each, are consistent with white noise: 

*r,s = 9r,s > -n r _< n, li.s.1 N s , (1) 

and N = 2n+l . 

Here X rjS refers to the data at day r in realization s; thus n = 15; N - 31; 
N s = 1 1’. The g r?s are random variables, uncorrelated with respect to both 
indices r and s, and satisfying 

E[g r>s ] = 0, for all r and s (2a) 

E[g r , s V,s'^ = f° 2 if r = r* and s = s’ 

10, otherwise . (2b) 

The notation E[ ] refers to the theoretical expectation or ensemble average 
over all possible independent realizations; the variance a 2 is a free parameter 
which can be chosen to fit the data. The precise meaning of consistency will 
be given in Eq. (6). 

The second hypothesis H] is that the data are consistent with a 
first-order autoregressive (Markov) process (Gilman et al . , 1963; Jenkins 

and Watts, 1968): 

V,s = a l^r-l,s + 9r,s > (3) 

where a^ and the variance (a 2 ) of g r «. are to be determined from a best fit to 
the data. 


The hypothesis of a second-order autoregressive process (H2) can be 
written as: 


x r,s = a l^r-l,s + a 2 x r-2,s + 9r,s 
where a^ , a 2 , and a 2 are to be determined. 


(4) 
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To better understand the nature of these low-order processes, it is 
helpful to consider the autocovariance function at lag k 

* £[X r>s » 

which measures the (average) influence of X at any day on its value k days 
later. The hypothesis of white noise implies that the value of X on any one 
day is completely uncorrelated with its value on any other day. Hence Yk * 0 
if k f 0. A first-order process describes "damped persistence," with Y|< 
smoothly decreasing from its value at k - 0 to a near-zero value some days 
later (see Fig. 3). It can be shown that second- and third-order processes 
can describe a damped periodic, as well as a purely damped, component so that 
Y(< resemb.les a sine wave with amplitude decreasing in time (see Fig. 4). 

The method of determining the adequacy of each of the hypotheses in 
describing the data is discussed in detail in Straus and Halem (1979). A few 
of the basic ideas are illustrated here using H? as an example. The "best" 
estimates of the parameters a^ and ag in Eq. (4) are obtained by minimizing 
the mean-squared error V, 

N s n 

( x r,s“ a l x r-l,s ' a 2 x r-2,s) 2 » ( 5 ) 

s=l r--n 


V(a-j , 


with respect to a] and a 2 » Let the estimates so obtained be written as a] and 
a 2 - Then, according to H 2 , the residual g r>s , given by 

X r> s “ ' a l x r-l 4 s “ a 2 x r-2,s = 9r,s » ( 6 ) 

should behave as white noise. In this manner, a test of all the hypotheses can 
be accomplished by use of a single statistical test for white noise. Such a 
test has been described in Straus and Halem (1979). A particular hypothesis is 
considered consistent with the original data if the appropriate residual is 
consistent with white noise. 

The lowest order autoregressive process that was consistent with the model data 
is shown in Fig. 1 (for ground temperature) and Fig. 2 (for sea level pressure). 
The hatching indicates stations for which none of the processes considered is 
consistent with the data. It is clear that white noise is not consistent with 
the model climatic noise, and that (except for the ground temperature at a few 
stations) this is also true of a first-order Markov process. There are large 
regions in the midwest and west for which a second-order process is consistent 
with both the sea level pressure and ground temperature. Note, however, that 
over portions of the northeast, none of the processes is consistent with the 
data for both variables. 

The almost total failure of the first-order process to capture the 
data is illustrated in Fig. 3, in which is shown the autocovariance function 
for sea level pressure calculated from the data and that implied by the "best 
fit" first-order process. For lags of 0 or 1 , the first-order process is in 


GROUND TEMPERATURE [ 31 DAYS ] 
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Fig. 1. 


Ground temperatures over the U.S. for an ensemble of eleven integrations of 31 days 
each. The lowest order autoregressive process that is consistent with the data is 
indicated. The hatching shows regions where white noise, first-, second-, and third- 
order processes are not consistent with the data. 






i pupi pressure over the U.S. for an ensemble of eleven integrations of 31 days 
each The lowest order autoregressive process that is consistent with the data is 
indicated? T^e hatching shows regions where white noise, first-, second-, and third- 
order processes are not consistent with the data. 
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Autocovari ance functions. The solid line is the autocovariance function obtained 

directly from the data for Austin, Texas. The dotted ^econfidellce 

function of the "best fit" first-order process. The dashed lines are confidenc 

limits on the first-order autocovariance runction and depend on the sample size. 



good agreement with the data. For larger lags, the sharp drop in the auto- 
covariance function and subsequent oscillation seen in the data are not 
present in the first-order process. The tendency for the observed auto- 
covariance to become negative quite rapidly is fairly general, and explains 
why the first-order process is not consistent with the data. Figs. 4 and 5 
show the autocovariance functions of the "best fit" second- and third-order 
processes. Since these processes can describe damped oscillations, they are 
much more consistent with the data. 

Two conclusions are suggested by these results. First, it may be 
misleading to employ simple stochastic model processes (e.g. , a Markov process) 
in explaining meteorological data (in this case climatic noise) because of 
the variation in stochastic properties with geographical region. Second, a 
possibility arises that for some regions none of the very simple univariate 
processes are consistent with observational or numerically simulated data, 
In these cases, it may be necessary to use more sophisticated models, in which 
the variable in question may depend on more of its own past history than the 
previous two observations (high-order autoregressive processes), on its value 
at other locations, or even on the values of other dynamically related varia- 
bles (multivariate processes). 
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Autocovariance functions. The solid line is the autocovariance function obtained 
directly from the data for Austin, Texas. The dotted line shows the autocovariance 
function of the "best fit" third-order process. The dashed lines are confidence 
limits on the third-order autocovariance function and depend on the sample size. 
(See Straus and Halem, 1979.) 



SIMPLE CLIMATE MODELS OF GLACIATION CYCLES 


(M. Ghil, K. Bhattacharya, C. Crafoord, E. Kalian, and H. LeTreut) 


An Energy-Balance Model with Time-Lagged Albedo 

The long-term energy balance of the climatic system can be expressed 
mathematically by a nonlinear partial differential equation (PDE) of parabolic 
type, with latitude x as the spatial variable and temperature T as the dependent 
variable. The crucial nonlinearity, appearing in the source term of the equation, 
is related to the physical dependence of global reflectivity, called albedo, on 
temperature, a* a (x,T). The function a(x,T) is assumed to be piecewise linear 
in T at constant x, leading to the existence of multiple equilibrium solutions 
(Ghil, 1976); some of these can be identified with possible climates of the 
earth. The delay t introduced (Ghil and Bhattacharya, 1979) represents the 
adjustment time of a to T, a *a (x,T(x,t-x)). 

The number of equilibria and their stability is studied as a function 
of the normalized average intensity of solar radiation, P . The behavior of 
time-dependent solutions is studied as a function of the delay t. Bifurcation 
to stable periodic solutions (Bhattacharya, 1979; Ghil and Bhattacharya, 1979) 
occurs with increasing t, followed by the appearance of multi-periodic and 
quasi-chaotic solutions (Bhattacharya, 1979; Bhattacharya and Ghil, 1978). This 
behavior seems to argue for the possibility of causes internal to the climatic 
system playing a role in glaciation cycles; no periodic external forcing is 
necessary in order to obtain stable oscillations of the system on the time scales 
involved. 

Energy-Balance Models with Ice-Sheet Dynamics and with Geodynamics 

In this work, the energy-balance equation for T is simplified by 
eliminating spatial dependence altogether. On the other hand, the complexity 
of the model is increased by adding an equation for the latitudinal extent 
of the ice sheets, yielding a system of two coupled, nonlinear autonomous 
ordinary differential equations (ODEs): c-jT * f(T,x)» C 2 & * g(T,£). This 
system has multiple critical poin f * 0 * g, one of which can be identified 
with the present climate of the earth. This point is a stable focus for a 
certain range of parameter values. As the parameter C 2 /c-|, corresponding 
to the delay in the simple partial differential equation (PDE) model, is 
increased, the focus becomes unstable, giving birth to a stable limit cycle 
(KallSn et al . , 1978). The amplitudes, period and phase difference between T(t) 
and i(t) for the periodic trajectories are of the order of magnitude given by 
paleoclimatological data on glaciation cycles (Katlen et al . , 1979). 

The model is being further refined to include the viscoelastic 
response of the lithosphere (crust) and asthenosphere (upper mantle) to the 
load of variable, size ice shee.ts (Ghil et al . , 1979). The system under study 
then becomes; qT - f(T,t), c £&= g'(T ,*>,"£), 03 i - h(t,t). Its critical points 
correspond to those of the preceding system, but their stability properties 
differ (Le Treut and Ghil, 1980). It is hoped that the oscillations of the 
latter model will have characteristics in even better agreement with the paleo- 
climatological record than the ones previously discussed. 
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SIMULATION OF THE ATMOSPHERIC GENERAL CIRCULATION ON MARS AND COMPARISONS 
WITH SPACE-PROBE OBSERVATIONS: A STUDY IN PLANETARY CLIMATOLOGY* 

(Y. Mintz) 


The atmospheric general circulation on Mars is simulated with a 
numerical model, which differs from the one used by Leovy and Mintz (1969, 
paper I) in that it takes into account the planet's large-scale surface 
topography. The new model also has somewhat greater horizontal and vertical 
resolutions, improved formulations for calculating the solar and infrared 
radiative heating, and a better way of parameterizing the vertical convective 
heat transfer within the planetary boundary layer. The simulation is for the 
time of southern hemisphere winter, which is the time of the Viking space- 
probe landings. 

In this simulation, the large-scale surface topography produces two 
kinds of forced standing waves: quasi -barotropic Rossby waves, which are 

kinematically forced by the topography upon the strong mid-latitude zonal flow 
of the winter hemisphere; and monsoon-like baroclinic circulations, which are 
thermally forced by strong heating over the summer hemisphere regions of high 
surface elevation. 

As in I, transient eastward moving waves are produced in the middle 
and high latitudes of the winter hemisphere. The dominant period and ampli- 
tude of these waves, as shown by the surface pressure and wind variations, 
agree with the Viking lander observations of the winter season transient waves: 
a period of about 3 days, root mean square (rms) surface pressure variation of 
about 0.15 mb, and rms surface meridional wind variation of about 5 m s"' 
(see Fig. 1). The dominant zonal wave number and phase speed of the transient 
waves also agree with what can be inferred from the Viking data: zonal wave 

numbers in the range 4 to 6, and phase speeds of the order 15 to 20 m s _ K 

The zonal ly and time-averaged fields of temperature and zonal wind, 
the mean meridional mass circulation, and the extent of the CO 2 surface ice 
cap, closely resemble those found in I, except for differences attributable to 
a southern hemisphere as against a northern hemisphere winter season. These 
features, therefore, do not depend on the presence or absence of surface 
topography, nor are they sensitive to the differences in the radiative and 
convective heating formulations of the two models. Both simulations show 
small horizontal temperature gradients and large vertical temperature gradi- 
ents from the summer pole to the tropics of the winter hemisphere, a narrow 
zone of large horizontal temperature gradient and small vertical temperature 
gradient in the mid-latitudes of the winter hemisphere, and a nearly three- 
dimensional ly isothermal state of the atmosphere over the winter polar cap. 
Both show a broad and intense thermally direct cell of mean meridional mass 
circulation centered on the winter side of the equator, and a narrow and weak 
thermally indirect cell on the polar side of the mid-latitude zone of large 
horizontal temperature gradient. In the winter hemisphere mid-latitudes, 
there is a very strong eastward (i.e., westerly) surface zonal wind, whose 

* Abstract and figures from paper by J. B. Pollack, C. B. Leovy, P. W. Greiman, 
and Y. Mintz, which is being submitted to J. Atmos. Sci. 
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Fiq. 1. Top panel : Simulated surface pressure, and simulated 

meridional wind component at tJs 0.89 (approximate height 
0.8 km), at 120W, 60S in southern winter. 

Bottom panel: Viking Lander 2 observed surface pressure, 

and observed meridional wind component at 1.6 m height, 
at 225W , 48N in northern winter. 
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maximum (15 m s -1 ) is south of the latitude of the strongest thermally 
indirect mean meridional circulation. The vertical shear of the zonal wind, 
as expected, has its maximum where the poleward temperature gradient is 
maximum and this, together with the surface zonal wind, determines the loca- 
tion and intensity of the upper level mean eastward jet stream (85 m s" 1 
at 42N). Within the tropics, in both hemispheres, the surface zonal winds 
are westward; decreasing in intensity with height on the winter side of the 
equator and increasing in intensity with height on the summer side of the. 
equator. Between the tropics and the pole of the summer hemisphere, there 
are several bands of alternating eastward and westward mean zonal winds, 
whose maximum magnitudes are about 10 m s -1 (see Fig. 2). 

The zonal ly and time-averaged simulated temperatures are in good 
agreement with the Mariner 9 remote infrared temperature soundings and with 
the Viking entry temperature measurements, except that the simulated lapse 
rates are larger than the observed ones. This may be because the model 
neglects the absorption of solar radiation by the dust in the Martian atmos- 
phere (see Fig. 3). 

The simulation shows a large rate of globally averaged surface 
pressure fall (0.011 mb sol" ) due to the condensation of atmospheric C0 2 
onto the surface in the winter polar region; and this agrees very closely 
with the Viking Lander surface pressure measurements. 

The simulated atmospheric heat transports are very large per unit 
atmospheric mass; but, because the mass per unit column of atmosphere is 
very small, the total atmospheric heat transports are small. Moreover, on 
the poleward side of the boundary of the condensing CO? polar cap, much of 
the poleward eddy heat transport is compensated by tne equatorward heat 
transport of the indirect mean meridional circulation, a circulation which 
is mechanically forced by the transient waves. Consequently, the total 
poleward heat transport by the atmosphere does not balance the radiative 
cooling of the winter polar region and, therefore, it does not keep the 
temperature of that region from falling to the CO 2 frost point. The polar 
cap of CO 2 ice that forms on the surface, in the simulation, reaches to 
latitude 39S. This is within about 1° of the boundary of the observed surface 
ice cap. 

The CO 2 surface condensation produces a net mass flux toward the 
condensing polar cap; and this produces a net eastward Coriolis torque, with 
its maximum at the edge of the condensing polar cap. It is this torque, 
together with the low-level eastward torque of the eddy-driven indirect mean 
meridional circulation, which determines the location and the great strength 
of the winter mid-latitude eastward surface winds. 
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Fig. 3. Meridional cross sections of the zonal ly and time- 
averaged temperature (°K). 

Top panel: Simulated temperature, in southern mid-winter. 

Bottom panel: Temperature derived from Mariner 9 infrared 

soundings, in southern late fall-early winter. 

The dashed lines and boxed-in numbers show the calculated 
geostrophic zonal wind (m/s) when the surface wind is 
neglected. 

The horizontal line, at p « 1 mb, indicates where the 
general circulation model has its upper boundary. 







ATLANTIC OCEAN SEA SURFACE TEMPERATURE ANOMALIES 

(J. Caiman) 


Since about three-fourths of the earth's surface is covered by 
water, any attempt to forecast the state of the atmosphere (or ocean) must 
take into account the physics of the air-sea interaction processes. The 
effect of this interaction is very strong in the upper hundred meters or so 
of the ocean where the temperature and depth of the "mixed layer" (so called 
because the temperature is constant from the sea surface to the bottom of 
the mixed layer) where the rapid fluctuations of sea surface temperature 
(SST) reflect the changes in the balance of turbulent fluxes of heat, mass, 
and momentum across the air-sea interface. 

In a recent study, Bunker (1976) reported measurements by ships of 
sea surface temperature, wind, humidity, and cloud cover were collected from 
the archives of the National Climatic Center. Using bulk transfer formulas 
for the parameterization of turbulent fluxes, Bunker studied the annual average 
(from 1948-1972) components of the heat budget of the Atlantic Ocean. Subse- 
quently, Bunker and Goldsmith (1979) have made available to other investigators 
their computed monthly values over the entire 25-year period of the components 
of the heat budget of the Atlantic Ocean. This data set has been obtained 
at the Laboratory for Atmospheric Sciences (GLAS) Modeling and Simulation 
Facility (GMSF), at NASA Goddard Space Flight Center, and is being used to 
study the behavior of SST anomalies in the Atlantic Ocean. It has been 
recognized that the usual ad hoc practice of specifying the mean monthly SST 
as a time-dependent boundary condition for an atmospheric general circulation 
model is an inadequate model of the air-sea interaction. In an effort to 
more correctly account for the physics of this interaction, ocean mixed-layer 
models, which describe changes in the mixed-layer properties (temperature, 
depth, velocity) as functions of the mass, momentum, heat, and turbulent 
energy budgets of the layer, have been developed. Several versions of the 
GMSF ocean mixed-layer model have been described elsewhere (Cane and Schopf, 
1978; Fung and Cane, 1978a, b). 

The 25 years of monthly measurements referred to above are being 
used to improve our understanding and modeling of the mixed layer in two 
ways. First, a descriptive study is underway to see how, on the basis of 
the measurements alone, the observed SST and wind stress, clouds, humidity, 
and calculated heating (from the bulk transfer calculations of Bunker) are 
statistically related to each other. Second, a study of the accuracy of the 
GMSF ocean mixed-layer model is being conducted. Previous comparisons of the 
predictions of the ocean mixed-layer model with observations have only been 
possible using the mean observed forcing and the mean observed response. It 
is of great interest for the purpose of constructing an interactive air-sea 
general circulation model to see how well the mixed-layer model predicts the 
observed anomalies (i.e., the departure from the mean) of SST given the 
observed anomalies in the forcing functions (heat, wind, stress). With the 
new data, such a test is possible and is now underway. 
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The annual and monthly means of SST, heating, and wind stress 
magnitude have been described (Bunker, 1976). The corresponding anomalies 
for January 1970 are shown in Figs. la-c. The time series of these anomalies 
at 35N, 65W are shown in Figs. 2a-c. It is surprising that correspondence 
between the three spatial patterns in Fig. 1 and between the three time varia- 
tions in Fig. 2 is not more apparent. If the heat absorbed by the ocean is 
unusually large and the wind stress is unusually weak, it is expected that the 
SST will increase. It can be deduced from Fig. 1, however, that this relation- 
ship is not generally true. Theory suggests that it is the rate of change of 
SST that responds to changes in heating, rather than the SST itself, and this 
relationship is being investigated. Regions of expected strong currents do 
not appear to explain the differences. Another physical argument suggests 
that sea. surface temperatures take some time to change and that, when the 
mixed layer deepens, it entrains cold water left behind the previous year. It 
is therefore expected that the ocean 'remembers* what happened at previous 
times. A look at the time series of Fig. 2 does not show an obvious time-lag 
correlation between the variables. The computed auto- and cross-correlation 
functions are shown in Fig. 3a for a point at 35N latitude, 65W longitude. 
The time scale of SST anomalies is several months, while those of heating and 
wind stress are less than 1 month and are unresolved in these data. (When 
linear trends are removed from the anomaly time series, the correlation 
functions approach zero at long lag times. ) The cross-correlation of SST and 
heating (Fig. 3b) shows that anomalous strong heating will lead to a warmer 
SST for several months, while a warm SST has no correlation with future 
heating. Similarly, the cross-correlation of wind stress with SST shows that 
strong winds will cause a cooler ocean for several months while a cool SST 
gives no knowledge of future winds. The cross-correlation of heating and 
wind shows no relationship between the two quantities except at zero lag where 
the high correlation is an artifact due to the use of wind in the bulk transfer 
formula used for computing the heating. Spatial variations of these correla- 
tions are being studied to identify special regions. An important test of the 
heat balance equation for the mixed layer is being made by computing the cor- 
relation of the rate of change of temperature with heating and wind stress. 
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Fig. 3b. Cross-correlations of SST, heating, and 
wind stress at 35N, 65W. 














ANALYSIS OF THE SURFACE HEAT BALANCE OVER THE WORLD OCEAN 

(S. Esbensen) 


During the period of the present NASA grant with Oregon State 
University, the major accomplishment has been the analysis of approximately 
25 years of ocean weather ship data at representative locations in the North 
Atlantic and North Pacific Oceans (see Fig. 1). A conclusion has been 
reached that monthly averaged surface meteorological data can provide heat 
flux estimates that are within 25 percent of the monthly mean heat flux 
estimates calculated from each observation. The weather ship data were 
obtained on magnetic tapes (in TDF-11 format) from the National Climatic 
Center, Asheville, North Carolina, and the Scripps Institution of Oceanog- 
raphy, La Jolla, California. 

To parameterize the heat fluxes, the bulk aerodynamic equations were 
used which are 


So - PCpC T V a (T s - T a ) (1) 

and 

LE 0 * pLC E V a (q s - q a ) (2) 


where S 0 and LE 0 are the surface sensible and latent heat fluxes, Cy and C E 
are the transfer coefficients, and V, T, and q are the wind speed, temperature, 
and -specific humidity, respectively. The subscripts "s" and "a" refer to sea 
level and near-surface (here assumed to be 10 m) values of atmospheric varia- 
bles, respectively. The parameters p, c p , and L are the near-surface air 
density, the specific heat of air at constant pressure, and the latent heat 
of condensation or evaporation. In this study, the transfer coefficients 
were determined from a modified version of the parameterization proposed by 
Benoit (1976) who used the suggestions of Deardorff (1972). 

Fig. 2 summarizes the heat flux computations for January and July. 
S 0 and LE 0 are monthly averaged values of the heat fluxes produced by sub- 
stituting 3-hourly observations into Equations (1) and (2), summing all of 
these flux estimates within a given month over the years of record, and 
finally dividing by the total number of observations. S 0 and IF 0 are the 
flux estimates obtained by substituting monthly averaged wind speeds, tem- 
peratures, and humidities into Equations (1) and (2). As shown in Fig. 2, 
the agreement between the two methods is better than 20 percent over a wide 
range of heat fluxes. Furthermore, this result may also hold for individual 
months. 
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. A comparison of long-term monthly mean heat fluxe 
(S 0 , LE 0 ) compared with heat fluxes computed from 
monthly mean data [Z 0 , LE 0 ) for (a) January and 
(b) July. 
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_ _ Jig. 3 s hows the comp arisons of V a (T s - T a ) with V a (Tc - T a ), and 
V a (q« - q^) with V a (q s - q a ) for twenty-three Januaries at snip D in the 
highly variable flux region over the North Atlantic Current. The overbars 
here refer to an arithmetic average of 3-hourly observations over the period 
of a month. Again, the agreement is remarkable. These results suggest that a 
combination of surface wind retrievals from satellite data and thermodynamic 
variables from ocean buoys may provide adequate data to monitor the monthly 
averaged surface heat fluxes in data sparse areas. 
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RESPONSE OF A GLOBAL MIXED LAYER MODEL TO GCM FORCING 
(I. Fung and M. Cane) 


This study reports the continued experimentation with an ocean 
mixed layer model before coupling it to the a nospheric GCM. 

The mixed layer model has been reported by Fung and Cane (1978). 
It has been shown to reasonably predict sea surface temperatures (SST) and 
mixed layer depths given climatological inputs of surface wind stress (t c ) 
and surface heat flux (Q c ). 

Climatological inputs are necessarily smooth, the value at each 
time step being linearly interpolated from means assumed at mid-month. 
Before the mixed layer model is coupled to the atmospheric GCM, it is 

desirable to study its response to high frequency forcing, an ideal source 
of which is the GCM. 

The GLAS GCM has been integrated for 4 months from winter initial 
conditions (Halem et al . , 1979). The ocean is noninteractive: SST is 

prescribed in time and space. From this GCM run, surface wind stress (t) 

and net surface heat flux (Q) were extracted every 6 h as forcing for the 

mixed layer model. Due to the lack of a boundary layer in the GCM, surface 
winds were linearly extrapolated from levels 8 and 9 to the surface. t 
was calculated from the surface winds using a drag law where the drag 
coefficient is a function of static stability. Q is the algebraic sum of 
the incident solar radiation, outgoing thermal radiation, sensible heat and 
latent heat of evaporation simulated by the GCM. Thus forced, the mixed 
layer model was integrated from climatological January initial conditions 
for 4 months, the duration of the GCM experiment. It is noteworthy that the 
mixed layer model remained stable throughout the 4 months to the high frequency 
atmospheric forcing. 

Fig. 1 presents the zonally averaged SST departure from climatology 
(aSST) for January. A positive (negative) ASST means that the simulated 
SST is higher (lower) than the climatological value. After 1 month, ASST 
is less than 1° over most of the ocean, with the largest errors at the 
latitudes of the ITCZ and the roaring 40' s. 

The trend of SST can be deduced easily from a comparison of the 
GCM-derived forcing with climatology. (See Fig. 2.) The much stronger 
winds at all latitudes act against the excess heat flux to produce very 
interesting results. 

North of 30N, even though Q>Q C (not enough cooling), t/t c ~3„5. 
Cold water entrained from below the mixed layer dominates over the reduced 
heat loss from the surface. ASST<0, but small; the errors nearly compensate. 
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Between ION and 30N, t/t c ~2.5., but heat is input into the ocean 
rather than lost from the ocean. ASST>0, but again the forcing errors 
nearly compensate. 

Between ION and 10S, Q>Q C (slight warming), t/t c *5. Heat is 
distributed into a deeper mixed layer. ASSTCO. 

Between 10S and 30S, Q>Q C (excess warming). Even though t/t c » 2.5, 
the strong winds cannot overcome the very stable stratification; a shallow 
summer mixed layer forms. ASST>0, but once again the errors nearly compensate. 

In the roaring 40' s, the initial conditions have a deeper mixed 
layer than 10S to 30S. The retreating mixed layer remains deep so that, 
despite Q>Q C , the excessive wind stirring results in cooling. This trend 
is continued through April where the SST departures of the order of 3°. 
(See Fig. 3. ) 

The ocean is noninteractive in the GCM experiment. With a GCM- 
mixed layer interactive experiment, an enhancement (reduction) of evapora- 
tion, sensible heat and thermal radiation loss from the sea surface where 
ASST is positive (negative) would be expected. The net heat flux into the 
oceans will be correspond!' ngtly reduced (increased), resulting in a reduc- 
tion of | ASST | . Thus, the SST departures obtained in this study may be 
viewed as an upper bound for the “errors" in SST expected in an interactive 
run. 


A further experiment has been performed where climatological wind 
stress and GCM heat fluxes were used as forcing for the mixed layer. The 
reduction in the wind stirring results in a shallower mixed layer making 
ASST more positive everywhere. 
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A ONE -DIMENSIONAL AIR-SEA INTERACTION MODEL WITH PREDICTED CLOUDINESS 


(D. A. Randall, J. Caiman, C-H. Moeng, and P. Schopf) 


Recent studies of the cloud-topped mixed layer have shown that 
stratocumulus clouds and the planetary boundary layer (PBL) both deepen as 
the sea surface temperature (SST) increases (for fixed divergence), but that 
for a sufficiently warm SST the cloud layer is abruptly destroyed as a result 
of penetration by the dry, warm overlying air. The influence of SST on PBL 
cloudiness is thus somewhat complicated. 

In most studies to date, the PBL structure has been computed for 
prescribed values of the SST and large-scale divergence. But, of course, the 
SST is partially determined by local solar insolation and long wave cooling, 
both of which are very sensitive to cloudiness. The surface turbulent fluxes 
of sensible and latent heat are also influenced by PBL cloudiness. 

To the extent that an increase in SST favors an increase in cloud 
depth, and to the extent that a deeper cloud favors a colder SST by reducing 
solar insolation, it can be anticipated that small fluctuations in SST ar,d 
cloud depth will be damped. But if a positive fluctuation of the SST is 
large enough to trigger the layer cloud instability, thereby reducing the 
cloudiness, then the SST fluctuation may be amplified. The coupled system 
thus appears to be stable with respect to small perturbations but unstable 
with respect to large perturbations. 

To study this and other aspects of local air-sea interactions, the 
pre-existing one-dimensional atmospheric cloud-topped mixed layer model of 
Randall (1980b) and the pre-existing one-dimensional oceanic mixed layer (OML) 
model of Fung and Cane (1978) and Cane and Schopf (1980) have been coupled. 
In order to provide for the radiative effects of clouds on the atmospheric 
turbulence and for the radiative warming and cooling of the sea surface, 
the model includes the radiation parameterization of Herman and Goody (1976), 
which includes the diurnal cycle of solar radiation. For the air, the prog- 
nostic variables of the model are the depth, wind, moist static energy, and 
mixing ratio of the mixed layer; for the ocean they are the depth, current, 
and temperature of the mixed layer. The model allows the possibility of a 
layer of nonprecipitating stratocumulus clouds, which promote convection and 
entrainment through phase changes and cloud-top radiative cooling (Randall, 
1979; 1980a, b). The surface fluxes of sensible heat, moisture, and momentum 
are determined through the bulk aerodynamic method. The free atmosphere and 
deep ocean soundings are prescribed as boundary conditions. Oceanic upwell- 
ing and atmospheric subsidence are varied parametrically. 

Fig. 1 shows results for the last 3 days of a 50-day test run of 
the coupled model. Diurnal oscillations of the SST, PBL depth, and OML 
depth are evident, although the SST and the PBL depth still show a linear 
trend with time, indicating that a longer run is needed to reach equilibrium. 
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The amplitude and phase of the SST and OML depth oscillations are reasonable, 
but the PBL depth oscillation has a smaller amplitude and different phase than 
anticipated on the basis of the results of Schubert (1976). This discrepancy 
has not yet been resolved. 

In the immediate future, shakedown runs of the model will continue, 
so that its sensitivities and overall credibility can be evaluated. Then the 
model will be applied to the study of cloudiness-SST feedback mechanisms. 
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Fig. 1. Sketch of meridional section of the steady x-independent 
response to easterlies and westerlies, looking eastward. 


A TWO-LAYER DYNAMIC AND THERMODYNAMIC MODEL OF THE UPPER OCEAN 

(M. Cane and P. Schopf) 


A new model has been constructed to simulate the thermodynamics 
and circulation of the upper ocean. The primary design goal is that the 
model be capable of state-of-the-art simulations of those processes that 
determine the sea surface temperature. The model has been kept simple so 
that it is inexpensive to run simulations over several years. The present 
fine-grid (As -24 km) model can be used in a stand-alone mode, either as a 
global model or to do experiments on a regional scale, or be coupled with 
the GLAS atmospheric GCM. 

The vertical structure of the model is depicted in Fig. 1. The 
upper layer, the mixed layer, is treated as a slab: Temperature and current 
velocities are taken as constant throughout. This is consistent with the 
bulk mixed-layer models now widespread (see Niiler and Kraus, 1977). The 
assumption that the layer is indeed well mixed is a good one for temperature 
(Fig. 2) but less so for currents. The depth of the mixed layer is variable 
and is determined by the thermodynamics , turbulent mixing, and divergence of 
fluid due to large-scale motions. 

The layer below the mixed layer is taken to extend to the permanent 
thermocline with the deep ocean approximated as being at rest (cf., e.g., 
Charney, 1955). In this second active layer, the currents are again assumed 
to be uniform with depth. This is the simplest structure that allows one to 
model the effects of static stability on mixed layer entrainment. In both 
active layers, the dynamics are governed by the primitive equations. The 
mixed-layer depth and temperature are also affected by local mixed-layer 
processes. These are modeled by a bulk formulation very similar to that of 
Niiler and Kraus (1977). Thus, in the model, as in nature, the surface 
temperature is influenced by local surface heating and wind stirring, hori- 
zontal advection, and upwelling. (Cooling due to upwelling is simulated by 
first having surface layer divergence thin the layer, thereby bringing the 
cold water below nearer to the surface where turbulent mixing processes 
entrain it into the mixed layer.) 

The model is being used to study the relationship between local 
thermodynamic and large-scale dynamical influences on sea surface temperature. 
The initial project focuses on determining the relative contributions to the 
seasonal heat budget of the Arabian Sea. Fig. 3 shows the sea surface tem- 
perature patterns predicted 25 days after the onset of the monsoon. This 
model run involves only the dynamically induced changes in surface tempera- 
ture, without any surface heat flux. Eddies have begun to form and move up 
the coast, causing strong gradients in the sea surface temperature. For 
comparison, the sea surface temperature on 6 July 1966, as determined from 
the NIMBUS II High Resolution Infrared Radiometer (HRIR), is shown in Fig. 4 
(Szekielda and LaViolette, 1970). It is apparent that the dynamically induced 
temperature anomalies are of the same scales and that the eddy signature 
observed from a satellite is similar to that produced by the model. 



Fig. 2. Typical variation of water properties with depth, November 1962. 

T = temperature, S = salinity, in parts per thousand by weight 
o - density anomaly, a = 1000 (density-1) where density is in 
g/cm^ and is measured at atmospheric pressure. (From H. L. Grant, 
A. Moilliet, and W. M. Vogel, 1968, some observations of the 
occurrence of turbulence in and above the thermocline, J. Fluid. 
Mech., 34* 443-448.) 
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Sea surface temperature observed from NIMBUS II satellite 
( HR I R ) during the SW-Monsoon of 1966. (From Szekielda and 
laViolette, 1970.) 








This initial study demonstrates the importance of nonlinear dynamics 
and high resolution for predicting even large-scale sea surface temperature. 
Further work with the model in its finer-grid version is planned to examine 
the low-latitude heat budgets and to prepare for coupling with the GLAS GCM. 
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RESONANT INSTABILITY AND MULTIPLE EQUILIBRIA IN OROGRAPHICALLY 
FORCED BAROCLINIC PLANETARY WAVE SYSTEMS 

J. G. Charney/MIT and D. M. Straus/National Research Council 


The general circulation of the atmosphere has often been regarded as 
consisting of a large-scale mean circulation, on which smaller scale 
transient waves are superimposed. The mean circulation, thought of as 
a Hadley circulation modified by topographic and thermal asymmetries of 
the lower boundary, is the single equilibrium of the atmosphere, albeit 
an unstable one. The smaller scale transients consist of waves and 
vortices which originate as instabilities of the mean flow and which 
interact with it. This description does not explain the presence of 
persistent, large-scale flow anomalies, such as blocking, which have so 
much to do with the atmosphere's observed variability. 

This presentation reports on research, both past and present, which 
suggests an explanation for the occurrence of persistent anomalies in 
terms of multiple equilibrium solutions. An equilibrium solution is 
a steady-state or periodic solution of the fundamental equations 
governing atmospheric dynamics. The existence of more than one 
solution is made possible by the asymmetries in the lower boundary, 
e.g.. mountains. If such solutions can be found which are stable to 
large-scale disturbances, one has an explanation for the existence of 
different large-scale flow patterns. If, in addition, these equilibria 
are unstable to smaller scale motions (as they will be in general), 
the transitions between the various large-scale patterns can also be 
explained. 

A simple model which incorporates these ideas is the spectral, quasi- 
geos trophic, barotropic, channel model studied by Charney and DeVore. 

The model is set on a mid-latitude beta plane, is severely truncated, 
but does include orography. Because of the truncation, there are only 
three degrees of freedom. One describes the zonal flow (which is forced 
directly), and the other two a single zonal wavenumber. Only one wave- 
length in the meridional direction ("first y-mode") is kept initially. 
The wave variables change due to advection by the mean flow, the beta- 
effect, friction, and an additional coupling to the zonal flow via the 
mountain (which has the same zonal wavenumber as the wave). The zonal 
flow is changed by forcing and dissipation, and also by interaction with 
the wave, again via the mountain. No Reynold's stresses are present . 

For small values of the forcing, only one steady state equilibrium 
exists, and it has non-zero wave components. It is also stable to 
perturbations of the first y-mode. However, for higher values of the 
forcing, three steady state equilibria exist. The two additional 
solutions are now stable, and the original one unstable. One of the 
stable states is a "high-index" circulation having a strong zonal flow 
and a relatively weak wave perturbation; the other is a "low-index" 
circulation having a weak zonal flow and a relatively strong wave 
perturbation. The unstable steady state has an intermediate flow. 
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This instability can be shown to be solely due to the interaction of the 
topographic wave with the zonal mean flow and is not a simple barotropic 
instability of the zonal flow. Numerical solutions of the three-component 
system, started near the unstable equilibrium, converge to either one of 
the stable equilibria, depending on exactly how the unstable equilibrium 
is perturbed. The presence and stability of these equilibria were con- 
firmed by numerical computations on a grid model having the same geometry, 
driving and dissipation as the spectral model, but having many more degrees 
of freedom. 

The interpretation can be made that blocking in the real atmosphere is an 
example of an equilibrium flow of the low-index type and that transition 
from it to the high-index state occurs via the instability of the intermediate 
state. 

The discovery of multiple equilibria and orographic instability in a baro- 
tropic system led us to investigate a corresponding baroclinic system. The 
important differences between the baroclinic and barotropic models include: 

(i) a description of the zonal flow and wave in each of two layers (for 
only the first y-mode) requires 6 variables; (ii) now it is the temperature 
gradient (or wind shear ) which is forced; and (iii) although Reynold’s 
stresses driving the average zonal flow in both layers (barotropic component) 
are again absent , the waves can directly change the zonal temperature gradient 
(baroclinic component) by causing heat fluxes to be present. 

As a consequence of (ii), the single equilibrium which exists for low values 
of the thermal driving is a simple Hadley solution, i.e. , a zonally symmetric 
solution with no lower level flow. (Hence, it is the same solution with or 
without mountains.) If orographic waves and a lower level zonal flow are to 
exist, they must be produced by the unstable growth of a small initial pertur- 
bation. Such an orographic instability (corresponding to the growth of a wave 
in place) is indeed found. However, as a consequence of (iii), it is not the 
only type of instability present; baroclinic instability in the usual sense 
(corresponding to traveling waves) Is also possible. The presence of two 
competing instabilities adds a great deal of richness to the behavior exhibited 
by the model. 

The simple 6-component model was studied in detail for wavenumber 3. For low 
values of the thermal driving, only the Hadley solution exists. The first 
stability threshold to be crossed as the driving is increased is that for 
orographic instability. This instability of the Hadley circulation leads to 
the existence of two stable, steady state equilibria with finite wave ampli- 
tudes ("wavy equilibria"). As in the barotropic case, the stable equilibria 
are both realizable in numerical integration. However, they both correspond 
to relatively high-index flows, i.e., ones with wave amplitudes small compared 
to the zonal flow. 

At higher values of the driving, the wavy equilibria themselves become 
baroclinically unstable, so that the system has three (including the Hadley 
solution) unstable steady state solutions. Numerical integration of the 
system leads to a periodic solution in which both the waves and the zonal 
flow vary in time. 
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One of the two wavy equilibria is* however, only very weakly unstable, and 
in this solution the wave amplitude increases dramatically with the external 
driving. At values of the forcing for which this solution begins to resemble 
a blocking (low-index) pattern, it becomes stable again. In this same regime, 
the Hadley solution becomes baroclinically (as well as orographieally) unstable, 
and the new instability leads to two new wavy equilibria, giving four in all. 

The two new equilibria correspond to very high-index flows, and are unstable. 

The blocking solution is stable, and the fourth wavy equilibrium unstable. 
Numerical solutions of the system, starting from either of the latter two wavy 
equilibria (the first two that appeared as the driving were increased from 
small values), converge to the stable blocking solution. However, starting 
from each of the high-index equilibria, the system converges to a distinct 
periodic solution. Thus, there are three stable solutions for fixed driving, 
one steady state and two periodic. * 

In both the barotropic and baroclinic systems, the stability of the various 
steady state solutions were also studied for a second meridional wavenumber 
(second y-mode), Many of the solutions that were stable to first y-mode 
perturbations became unstable to the second y-mode, and Integration of the 
appropriate equations with both y-modes led to the discovery of new periodic 
and also aperiodic solutions. Again, multiple solutions were found for 
fixed values of the driving. 

* 

These findings are regarded as useful primarily for their heuristic character, 
not for their detailed explanations of specific phenomena. They suggest that 
the study of steady state and periodic solutions, as well as consideration of 
the instabilities of these solutions, will become a useful tool for investiga- 
tion of large-scale atmospheric phenomena. We believe that this approach can 
lead to a better understanding of large-scale variability, predictability, and 
climate. To know that a system is close to a stable or metastable equilibrium 
is to know that it will remain so for a time. To know that the system is in a 
state of transition is to know that it will change more rapidly and be less 
predictable. Climate itself becomes a question of distributions among possible 
equilibrium states, and climate variation a matter of how boundary changes lead 
to altered distributions. 
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June 28, 1979 


THE EVOLUTION OF THE EXTRATROPICAL CYCLONE AND 
ITS INTERACTION WITH LARGER AND SMALLER SCALES 

D. Johnson/The University of Wisconsin 


The forcing of the mass circulation and angular momentum transport for two 
district types of cyclone development is summarized for the leeside 
development of the Alberta Cyclone of March 30-April 2, 1971. The first 
type of development during the initial 48 hours of the cyclone's lifetime 
was characterized by an inward branch of a mass circulation in the upper 
isentropic layers and an outward branch in lower isentropic layers. The 
upward branch of the mass circulation was forced by the combination of 
negative pressure and Inertial torques while the outward branch was forced 
by the combination of positive pressure and inertial torques. The transfer 
of storm absolute angular momentum was downward through the action of the 
torques with the result that the "spin-up” of the relative circulation of 
the extratropical cyclone occurs in the lower troposphere. Figure 1 from 
Johnson, Wash and Petersen (1976) shows that the tendency of specific storm 
angular momentum is positive in the lower layers, while it is negative in 
upper layers. During this type of development, precipitation is light and 
the flow is nearly adiabatic. Through a decomposition of the lateral 
transport processes into geostrophic and ageostrophic modes, the results 
indicate that development occurs primarily through the modes identified 
with baroclinic amplification of waves and that the cyclone develops 
in the low troposphere through larger scale forcing in the presence of 
topography. 

The second type of development occurring during the laBt 48 hours of the 
cyclone's lifetime was characterized by an inward mass transport in the 
lower isentropic layers and outward transport in the upper isentropic 
layers. The inward branch of the mass circulation was forced by an excess 
of the negative pressure and frictional torques over the positive inertial 
torque while the outward branch was forced by the excess of the positive 
pressure torque and eddy angular momentum convergence over the negative 
Inertial torque. During the transition stage between the two types of 
development, the precipitation within the cyclone's circulation Intensifies, 
and angular momentum is transferred upwards in association with ascending 
dlabatic trajectories and the release of latent heat along with the action 
of pressure torques. The qualitative contrast of the areal extent of the 
precipitation processes during the lifetime of the cyclone is portrayed in 
the three satellite pictures on three successive days (Figure 2), The first 
of the three geosynchronous images is taken during the second day of the 
cyclone's lifetime and displays little convection, while the last two taken 
during the third and fourth day show the extensive convection that developed 
within the cyclone circulation. The excess of the inward transport of 
angular momentum towards the cyclone's axis of rotation over the outward 
transport and the transfer of angular momentum upwards serves to "spin-up" 
the relative circulation at all levels and rapid occlusion follows. Note in 
Figure 1 that the tendency of specific storm angular momentum is positive at 
all levels. The partitioning of the lateral transport processes reveals 
that this type of cyclone development primarily occurs through the ageostrophic 
component of the mass circulation and that the mass circulation is forced 
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Figure 1. Tendency of specific storm 
angular momentum for Alberta Cyclone 
of 30 March - 2 April 1971. 








by physical processes identified either with the cyclone vortex or smaller 
scales embedded within the cyclone. This period of growth of the small- 
scale cyclone from a secondary vortex embedded within an amplifying baro- 
clinic wave to a large-scale occluded vortex within a large-amplitude wave 
is one where small-scale processes associated with the boundary layer and 
latent heat release through convection and frontal circulations force 
larger scale development. 
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STRATIFIED FLOW OVER (AND AROUND?) ISOLATED TOPOGRAPHY 
R. B. Smith/Yale University 


July 5, 1979 


The theoretical study of the airflow over mountains began with Lyra (1943) 
and Queney (1947). Since then, the subject has received a great deal of 
attention, in part, because the mathematical problem is well posed and, in 
part, because the flow can be considered steady. Still, in spite of these 
simplifications, the mountain flow problem has the basic elements of all 
problems in dynamical meteorology, namely the importance of the degree to 
which the flow is hydrostatic and/or geostrophic. The equations used by 
Lyra and Queney illustrate this point (also see Smith, 1979a). 


,.3u" . , 3P' 

pV 8x ~ P fv " “ * a ~ 


,.3v' _ . 3P * 

+ pfu “ - w 


pU, 


3w/ 

3x 


3P' 

3z 


BP 


3u' 3v' , 3w' n 

**“■■*• ***** m n 

3x 3y 3z 


U 3p ' 1 xt 2 , 

P Bx g 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


Noting the importance of the horizontal scale L on the advective terms in 
1, 2, and 3, we can reason (with benefit of hindsight) as follows. For 
flow over narrow mountains, the Vertical force balance will involve only 
the first two terms in (3) but, as L increases, the buoyancy force becomes 
relatively more important and eventually the terms on the right-hand side 
of (3) will dominate and hydrostatic balance is achieved. 

As the horizontal scale continues to increase, equations (1) and (2) 
change from a balance between advective acceleration and pressure gradient 
force to a balance between Coriolis force and pressure gradient force 
i.e., geostrophic balance. * 
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Generally speaking, as the buoyancy forces increase in importance and the 
flow becomes more nearly hydrostatic, the fluid becomes able to support 
internal gravity waves and the disturbance produced by the mountain can 
extend to great heights (Figure 1). [Several GOES film loops are shown 
indicating the high level cirrus formed in vertically propagating mountain 
waves over the Appalachians and the Front Range.] The theory of mountain 
waves also predicts the existence of a pressure difference across a 
mountain caused by mountain wave generation. Recently, this pressure 
difference has been directly measured using sensitive microbarographs 
placed on each side of a narrow mountain ridge (Smith, 1978). 

Most of the theoretical work on mountain flows has been restricted to 
two-dimensional flow over an infinite ridge but recently solutions have 
been obtained for flow past an isolated mountain (Smith, 1980). These 
solutions show, for the first time, the tendency for the low-level flow 
to be deflected to each side of the mountain and for the laterally dis- 
placed air to be replaced by the descent of potentially warmer air from 
aloft. This pattern of surface flow is associated with the generation 
of a 3-D wave pattern aloft. 

As we proceed to larger (i.e., synoptic) scales of motion, the Coriolis 
force becomes important. Eventually, when the Rossby number (U/fL) is 
small, the Coriolis force will dominate, quasi-geos trophic theory becomes 
valid, and the ability of the fluid to support internal gravity waves is 
destroyed by the Coriolis resistance to horizontal divergence. Recent 
theoretical and observational work (Smith, 1979b) has shown, however, 
that even for mountain ranges as wide as 300 km or so, mountain waves 
continue to play an important role (see Figure 2). It seems then that 
the quasi-geos trophic approximation should not be used in the analysis 
of flow over moderate-scale orography as the mountain waves are thereby 
lost; and with them go several important qualitative features — the 
propagation of the disturbance to great heights, the pressure difference 
across the mountain, and preference for flow to go around, rather than 
over, the mountain. 
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Figure 2. Three examples of the perturbation to the Surface wind 
and pressure fields caused by orography (a) The Scandinavian 
mountains , March 24,1953 0900z (from Spinnangr and Johansen } 1954), 
(b) Iceland ( from Bergfhorsson and Sigurdsson, 1976) and (c) Madagascar, 
July 18, 1964- (from Donque, 1975). The length scale of these flows 
( a few hundred km) is large enough so that the Coriolis force is 
important but the pressure asymmetry- high pressure upstream and 
low pressure downstream - seen in the figure ts evidence that the 
Coriolis force is not so dominant that it prevents the generation 
of mountain waves. 
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. A VARIATIONAL TECHNIQUE FOR THE ASSIMILATION OF METEOROLOGICAL DATA 

[ J. R. Bates/Irish Meteorological Service 


A variational method is proposed for global assimilation of meteorological data 
in which information from the previous forecast is used more fully than in 
existing methods. In most existing data assimilation schemes, the previous 
forecast is used only in providing a first guess for the analysis, while the 
initialization takes no further account of forecast information. Included in 
that category is the method based on nonlinear normal mode initialization. 

In the scheme proposed here, the forecast tendencies of vorticity and divergence, 
together with the forecast vertical motion and diabatic heating (all smoothed 
over the final three hours of the forecast), are used to modify the analyzed 
mass and wind fields, with the aim of preserving the maximum amount of useful 
historical information and of ensuring a smooth evolution of the computed 
atmospheric flow. 

Let (i|>°, <j>°) be the analyzed stream function and geopotential and let x° be the 
velocity potential derived (at the first iteration) from the forecast vertical 
motion. Let (°^, Otj, &j) be precision moduli associated with ($°, X°, $°). 

We minimize the functional 
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over the globe, where (>p, X, <{>) are functions to be determined. As constraints 
on the functional (1), we choose forms of the vorticity and divergence equations 
in which the tendencies and the nonlinear terms are prescribed: 
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The superscripts "F" denote foi^cast values. The nonlinear terms in a^ and &2 
are derived from the observed (analyzed) values at the first iteration, but 
are updated using the new values (ij>, X. 40 at succeeding Iterations. 

After the first iteration X° is derived from a vertical motion given by the 
thermodynamic equation 
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where T is derived from the hydrostatic equation 
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In this way, all the equations of motion enter the iterative procedure. 

To minimize the functional (1) subject to the constraints (2) and (3), We define 
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where A^, A2 are undetermined Lagrange multipliers. The Euler-Lagrange equations 
for this system are found to be 
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with similar equations where X, i replace V». 


Applied to (5), these equations become 


2o i v2 «'-r) + Jj- (AjB) + Lap - o 


( 6 ) 


2ct. 


v <X~X°) +8. (A 2 B) - Lap- 0 


(7) 


2a 3 (<j>-0°) + v 2 x 


( 8 ) 


where we have assumed (a^, (*2> a 3) are constants and 
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The equations (2), ( 3 ), (6), (7 ), (8) in (ip , X, A A, , A„) can 
expanding in spherical harmonics. Letting ^ ^ 


be solved by 
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where p*sin0 and p (p) are the normalized Lagrange functions, we find, 
after eliminating from (3) by means of (8) and using some standard 
recurrence relationships, 
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These equations can be written 


AX * B 


where A is a real square matrix and B is a complex column vector of order 4(J+1). 
The solution X, when iterated over all levels until convergence occurs, gives the 

adjusted motion field which is consistent with the tendencies of the previous 
rors^tis n • 
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July 19, 1979 


ON THE GENERATION OF LARGE-SCALE DIVERGENT WINDS AND RELATED ENERGETICS 

J. Paegle/The University of Utah 


Burger (1958) demonstrated that the ultra-long waves of the atmosphere 
should obey a simplified vorticity equation of the form 


f M - 

ap 


where f and 8 are the Coriolis parameter and its meridional gradient, 
respectively, co is the vertical motion and U is the meridional wind. 

Integration of (1) over pressure gives 

P 

»(P) - J M dP ^ 


In mid-latitudes, the integrand of (2) has a magnitude on the order of 
10~ 2 day -1 , for u on the order of 10 m/s. Thus, for essentially barotropic 
ultra-long waves, (2) suggests 


to (surface) * 0 (100 mb/day) 


while for ultra-long waves for which u changes sign with height 


co (surface) << 100 mb/day 


Therefore, for the barotropic case, parcels of air that follow ultra-long 
wave trajectories on the surface of the earth would have pressure variations 
that greatly exceed observed values, except in the vicinity of high mountains. 
It may be concluded that barotropic planetary scale waves are possible only 
in the presence of high mountains that have a significant projection upon 
ultra-long wave components. 

Figure 1 presents 200 mb and 850 mb meridional flows summed over longitudinal 
wavenumbers 0-3 at various latitudes. The data is from the DST period 
February 1, 1976-March 5, 1976, and is described in more detail by Paegle, 
et al. (1979). The fundamental conclusion is that from about 28N, southward 
into most of the Southern Hemisphere, the ultra-long waves are more similar 
to vertically reversing monsoonal circulations than to nearly barotropic wind 
systems shaped by topography. In the mid-latitudes of the Northern Hemisphere, 
the present results support the important role of topography found in many 
other studies. 
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Figure /, Long wave contribution to the meridional flow component 
at indicated latitudes at 200 mb ( solid curves)) and 850 mb ( dashed 
curves) for the time average of the winter 1976 DST data. The 
correlation coefficients between the two levels are indicated by 
r values. 
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Figure 2 displays the 200 mb divergent meridional wind for the same data set. 
It shows strong upper level outflow from the heavily precipitating regions 
over the Amazon basin, equatorial Africa, and the tropical Western Pacific, 
there is some similarity of these areas to the 200 mb ageostrophic meridional 
flow from high to low pressure shown in Figure 3. These poleward divergent 
and ageostrophic flows appear to emanate from strong convective heating zones 
calculated by Herman and Schubert (1979) for this data set and are associated 
with accelerations of the subtropical jet stream displayed in Figure 4. 

The meridional energy flux of the time-averaged fields is shown in Figure 5. 
It indicates poleward energy fluxes from about 40S to about 20N, at 200 mb, 
Indicating that the stationary features of the subtropical jet stream may be 
influenced by the distribution of tropical latent heating. 

These results suggest that longitudinal gradients of tropical heating play 
a more direct role in the forcing of long tropospheric, waves than is usually 
indicated. Further study should be done with forthcoming FGGE data analyses 
to test the generality of the conclusions. This could be particularly 
important if the short-term predictability of energy fluxes and divergent and 
ageostrophic flows is as sensitive to tropical, and even cross-equatorial, 
effects as suggested in modeling studies presented by Paegle, et al. (1979) 
and by Lewis (1979). 
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entire 2. Divergent meridional ;r .id component at 200 mb, averaged over 
the winter ,976 OS' period. Shaded regions indicate southerl ies 2 m/s 
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Figure 3. Cross-isoba^ic meridional wind component at 200 mb, 
averaged wer the vinter iq?6 DST period. Shading indicates regions 
where the ^poleward ageostrophic values exceed A m/s. The region 
between 6 N ond Q S is not analyzed because o-f the poor !y Jetermir ied 
qeosrtrophic wind at these latitudes. 
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. 5 Meridional energy flux as estimate from the correlation 
potential and meridional flow far the averaje 200 mb DST data. 



July 26, 1979 


WIND WAVE DEVELOPMENT 
E« Mollo-Christensen/MIT 


A presentation Is made of observations of wind waves on the ocean and as 
observed In the laboratory during a recent experiment performed jointly with 
A. Ramaonjiarisoa at 1' Institute Mecanique Statistique de la Turbulence at 
the University of Alx-Marseille 11. 

The new information shows the existence of subharmonic transitions in the wave 
field, as predicted by the analysis of Longue t-Higgins (1978). This subharmonic 
instability is followed by crest-pairing, a process by which one crest simply 
overtakes the previous crest and disappears (Figure 1). Through this process, 
the wave loses momentum, and this momentum has to be transferred either to the 
mean flow as drift current or to waves of larger scale. If the momentum transfer 
takes place at a favorable phase of the larger scale waves, the larger scale 
waves will grow by absorbing momentum and thus also some wave energy. 

This process is demonstrated in experiments on group collision in the absence of 
wind, in experiments with group collision in the presence of wind and naturally 
occurring wind waves, and also by observation of naturally occurring field. 

The process of wind group formation is also discussed. First, group formation by 
crest-pairing, then group merger because group of different length formed through 
crest-pairing and merger tend to have slightly different group velocity. The 
group lengths, measured in terms of wavelengths, tend to cluster around the 
Fibonacci numbers (Figure 2). When the groups become very long, through merger, 
they tend to split again; the result is a recurrence process, similar to the 
Fermi-Pasta-Ulam recurrence. 

In conclusion, while nonlinear waves tend to propagate and form groups according 
to the predictions of third -order theory, including following the nonlinear 
Schroedlnger equation for modulations, the presence of crest-pairing violates 
the conservation of crests criterion, so that, to the modulation equations one 
will have to add transition rules, yet to be defined. It also appears that 
internal wave packets, as observed on the seasonal thermocline, will tend to 
cluster in lengths around the Fibonacci multiples of the wavelength, as suggested 
by satellite observation. Crest-pairing may also be a process of importance in 
other kinds of wave fields, possibly including topographic Rossby waves. 

Yet another observation is that the sudden occurrence of subharmonic fluctuations 
also takes place in the air boundary layer over the wave field. This could well 
be caused by vortex pairing, as observed in free shear layers by Winant and 
Browand (1974). While the results are preliminary, they suggest that strongly 
nonlinear processes that occur very rapidly, and thus will tend to appear 
negligible in a scaling analysis, may still be of crucial importance. 
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August 9, 1979 


NUMERICAL MODELING OF STRATOCUMULUS IN THE BOUNDARY LAYER 
J. W. Deardorff /Oregon State University 


Much of the motivation for studying the turbulence in a boundary layer which 
contains stratocumulus in its upper portion lies in predicting the entrainment 
rate. The mixed-layer growth rate is the difference between the entrainment 
rate and the large-scale vertical velocity at the mixed-layer top. The 
boundary layer height is then an integral of the growth rate over time. 

The stratocumulus-capped marine layer often situated just off the west coast 
of the US, for example, may, because of subsidence, frequently be too shallow 
to cross the coast range of western Oregon in a sea-breeze situation in summer. 
The interior valleys may swelter in a heat wave until the marine layer can 
grow another 100pm deeper and flow through the passes in the coast range. 
Satellite photographs clearly illustrate the problem. 

The three-dimensional numerical model of the clear boundary layer developed 
by Deardorff ( Boundary-Layer Meteor. , 1974) was modified and applied to a 
boundary layer that could be saturated in its upper portion. Recently, the 
results have been analyzed for the case of the "dry cloud" with or without 
strong surface heating but with radiative cooling at its top, for a real 
cloud occupying the upper half of the boundary layer with or without 
radiative cooling at the top, and for a cloud-filled boundary layer with 
radiative cooling at its top. 

In the dry-cloud cases, the cloudtop cooling was inserted at a fixed height 
within the upper mixed layer but below the mean capping inversion or stable 
layer. The enhancement of buoyancy flux in the upper mixed layer was very 
noticeable, and it was clearly seen how the mechanism alone, in the absence 
of any surface heat flux, could generate the turbulence that would keep a 
mixed layer well mixed. 

In the numerical simulations of real stratocumulus, the influences of cloud 
water were very important in enhancing the buoyancy flux in the cloud layer. 
This was due to two reasons: (1) the upward moisture flux in the cloud 

layer caused a positive buoyancy flux there because of the Clausius-Clapeyron 
relationship, and (2) evaporation of cloud droplets into air entraining down 
into the cloud top caused this entrainment to proceed more quickly and easily. 
Thus, the negative buoyancy flux of entrainment that was observed to occur in 
the cloudtop region rapidly switched over to a positive buoyancy flux at 
slightly smaller heights. A strong positive moisture flux was also observed 
in this region, due to the entrainment of dry air. Consequently, the mean 
moisture content of the mixed layer tended to decrease with time in the deep 
mixed layers studied, and cloud base tended to rise with time, even when the 
evaporation rate from the surface was appreciable. 

In the real cloud simulations, the addition of radiative flux divergence 
(longwave) to the upper 50ym of cloud within each grid column produced no 
appreciable change in shape of the buoyancy-flux profile and turbulence 
energy. This was apparently due to the fact that the cloudtop hummocks were 
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located within the horizontal mean capping inversion, and the radiative 
extinction scale utilized (less than 50pm) was smaller than the depth of this 
mean capping inversion. However, the opposite result probably would have been 
obtained if the cases studied had been for shallow mixed layers and if the 
radiative extinction length scale had been considerably larger (100-150pm). 

With or without radiative flux divergence at cloud top, the entrainment rates 
obtained from the numerical simulations with cloud-water present, even after 
normalization by the mixed-layer convective velocity scale, were found to be 
about an order of magnitude greater than for, a clear mixed layer with comparable 
capping-inversion strength . This result is ascribed to the powerful influence 
of evaporative cooling acting on air which entrains down into the mixed layer. 

Unfortunately, due to model truncation error at and above simulated cloud top, 
many of the conclusions on the entrainment rate and entrainment buoyancy flux 
were mote qualitative than quantitative. Future numerical simulations of 
this nature could benefit substantially by using a much finer vertical grid 
than 50ym, at least in the vicinity of the mixed-layer top. 


August 16, 1979 


CLIMATIC EFFECT OF AN INCREASE OF C0 2 IN THE ATMOSPHERE 
S. Manabe/Geophysioal Fluid Dynamics Laboratory, NOAA 


A recent study of the climatic impact of future anthropogenic increase in 
the atmospheric CQo-concentration is described. For this study, a global 
model of climate with realistic geography and with seasonal variation of 
insolation was constructed. The model consists of: (1) a spectral model 

of the atmospheric general circulation, (2) a simple mixed-layer ocean 
model with a uniform thickness, and (3) a heat- and water-balance model 
of the continents. It is encouraging that this climate model successfully 
reproduces many of the basic characteristics of the seasonal variation in 
the geographical distribution of atmospheric temperature. The climatic 
effects of a (^“increase are evaluated by comparing two model climates, 
with the normal and four-times-the-normal concentration of carbon dioxide. 

In conclusion, the warming of the model atmosphere in response to a 
quadrupling of CC^-concentration Is most pronounced in the lower tropo- 
sphere of high latitudes due to the poleward retreat of the highly 
reflective snow cover and sea ice. The warming over the Arctic Ocean 
and the surrounding regions is at a maximum in early winter and is at a 
minimum in summer. Thus, the amplitude of the seasonal temperature 
variation reduces markedly. The reduction of sea ice thickness is 
responsible for the aforementioned seasonal asymmetries in warming. 
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August 23, 


SOME RESULTS WITH THE UCLA GENERAL CIRCULATION MODEL 
M. Suarez/UCLA 


The latest version of the UCLA general circulation model (GCM) was discussed. 
Simulations made with earlier versions of the model showed that, although 
the model compares favorably with the state of the art, it had several 
deficiencies. Some deficiencies were particular to the UCLA GCM; others 
were common to most GCM's. In an effort to improve the model, both for 
general circulation studies and for numerical weather prediction, two major 
design changes have been made: 

(1) A potential ens trophy-conserving advection scheme in the 
equation of motion has replaced the previous scheme, which only conserved 
enstrophy for non-divergent flow. It is hoped that this scheme will improve 
the simulation of topographically forced motions. 

(2) The treatment of the planetary boundary layer (PBL) has been 
dramatically modified by making the predicted PBL-top a coordinate surface. 
These and other more minor changes were discussed. 


August 30, 1979 


ON TROPICAL CIRCULATION AND CLIMATE ANOMALIES 
S. Hastenrath/The University of Wisconsin 


Variations of circulation and climate have been studied for a portion of 
the tropical belt extending from the Americas across the Atlantic to Africa 
and the Indian Ocean. Ship observations during 1911-1972 compiled as individual 
monthly means for one and five degree square areas provided the major data 
basis. This ensemble of observations is supplemented by conventional surface 
climatological and upper-air records at land stations and series of lake levels 
and river discharge. 

A background climatology has been established for reference and has included 
the publication of atlases of the climate and circulation, and the oceanic heat 
budget, of the tropical Atlantic and Eastern Pacific (Hastenrath, 1977a, b,c,d; 
1979; Lamb, 1977; Hastenrath and Lamb, 1977a, b; 1978a, b; Hastenrath and Guetter, 
1978). A two-volume climatic atlas of the Indian Ocean is in press (Hastenrath 
and Lamb, 1979). 

Mechanisms of climatic hazards have been studied with emphasis on the Central 
American-Carlbbean region, the Ecuador-Peru coast, Northeast Brazil, Subsaharan 
Africa, and the Indian Subcontinent (Hastenrath, 1976; 1978; Hastenrath and 
Heller, 1977; Covey and Hastenrath, 1978; Kaczmarczyk and Hastenrath, 1979; 

Lamb, 1978a, b; Hastenrath and Wendland, 1979). Even during extremely deficient 
or abundant rainy seasons in a given region, departures in the large-scale 
circulation are small, albeit distinct. Statistical significance of departures 
can be tested against the climatic variability for the 60-year base period. 

It has further become apparent that patterns must be viewed in context over 
very large areas. With this provision, characteristic departure patterns can 
be identified in the large-scale circulation. 

For example, during deficient rainy seasons in the Central American-Caribbean 
region (Hastenrath, 1976), the equatorward flank of the North Atlantic high is 
strengthened, the Northeast trade winds speed up, albeit in a band further 
South, and marked negative sea surface temperature anomalies extend from the 
realm of the cold Canary current across the North Atlantic in a broad zonal 
band between 10-20N; a tendency for opposite sea surface temperature departures 
is observed over much of the tropical South Atlantic; and the Eastern Pacific 
is anomalously warm in similarity to El Nino events which have a tendency to 
coincide with Central American-Caribbean droughts. Stratification with respect 
to droughts in Subsaharan Africa yields similar departure configurations (Lamb, 
1978a, b). Abundant rainy seasons are concomitant with approximately inverse 
circulation anomaly patterns. 

Rainfall in Northeast Brazil (Hastenrath and Heller, 1977; Covey and Hastenrath, 
1978) is related to the southernmost seasonal migration of the near-equatorial 
confluence axis and convergence band over the Western Atlantic. During drought 
years, these quasi-permanent circulation features stay comparatively far North, 
concomitant with an equatorward expansion of the South Atlantic and a poleward 
retraction of the North Atlantic subtropical highs. This is accompanied by 


positive sea surface temperature departures over much of the tropical North 
Atlantic and negative anomalies in the South Atlantic and the adjacent equatorial 
waters. Again, abundant rainy seasons in Northeast Brazil are characterized by 
broadly inverse atmospheric-hydrospherlc departure patterns. 

A signflcant negative correlation between Northeast Brazil rainfall and Eastern 
Pacific sea surface temperature (Covey and Hastenrath, 1978), as well as an 
apparent linkage between events on the Indian Subcontinent and the Americas, 
seems to be related to a large-scale mass exchange of the Southern Oscillation 
type, but requires further study. On the whole, there is a strong spatial 
correlation in annual rainfall totals in the zonal direction over thousands of 
kilometers whereas, in the meridional direction, correlations rapidly fade out 
and reverse sign. 

Map analyses, based on stratification with respect to extreme hydrometeorological 
events in key regions, are complemented by principal component analysis and 
spectral techniques (Hastenrath, 1978; llastenrath and Wendland, 1979; 

Kaczmarczyk and Hastenrath, 1979). Principal component analysis (Hastenrath, 
1978) identifies preferred modes of departure patterns in the pressure and 
sea surface temperature fields. These corroborate the various regional strati- 
fication studies, in that similar spatial patterns are obtained and the time 
series of principal components are significantly correlated with hydrometeoro- 
logical index series in key regions. Spectral techniques (Hastenrath and 
Wendland, 1979; Kaczmarczyk and Hastenrath, 1979) were used in an attempt at a 
preliminary inventory of the preferred time scales of variability in the regional 
climates and of the large-scale spatial linkages of regional extreme events. 

Major departure characteristics of anomalous rainy seasons are already apparent 
at the height of the preceding winter circulation, the inertia of the atmosphere- 
ocean system, thus offering the potential of seasonal foreshadowing. Our work 
so far has been essentially diagnostic. Prospects for prognosis have emerged 
particularly from departure fields and time series of pressure, wind, and 
sea surface temperature in sensitive ocean areas. Timely availability of 
predictors will he crucial in operational prediction schemes. Remote sensing 
by satellites could become highly useful in such endeavors, as satellite data 
banks geared at the needs of climate diagnostics and prognosis come into being. 
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A SIMPLE MECHANISM FOR BLOCKING 

E. Ka lnay-Ri va s / NASA GSFC and 
L."*0. Merkine/Technion-Israel Institute of Technology 


In this study, air *,how that atmospheric blocking can occur as a result of 
local resonant i,nleractlon between forced Rossby lee waves. In order to 
study the simples; possible model of local generation of blocking, we use the 
barotroplc vorticlty equation on a P-plane open channel. Unperturbed zonal 
flow entering the channel infringes on an isolated mountain (Figure 1)* and 
produces steady state Rossby lee waves (Figure 3a). We then allow interaction 
between this basic solution and Rossby waves produced by another source of 
localized forcing. 

In the first part of the paper, we find by means of an analytical study of 
the far field solution, that friction plays a singular role in the dynamics 
of Rossby waves. In the presence of friction, Rossby lee waves generate a 
rectified current downstream of the orography , which does not vanish if we 
let friction go to zero. The generation of Rossby lee waves, and the 
orientation of the rectified current depends on the width of the channel Ly, 
through the r^tio between the minimum meridional wave number (1 $/ l^) and 
$/U. Figure 2 shows a schematic representation of the steady state solution 
as determined from the far field analysis. The characteristics of the near 
field (dashed in Figure 2) are inferred from the perturbation energy equation. 

Multiplying the potential vorticlty equation 


(f^+Ul+^V) (vV + Py + f D jH). - ^ * a> 

by and integrating over the open channel, a perturbation kinetic energy 
is obtained: 


9KE’ 

~5t~ 


D' + B' 


( 2 ) 


where G' , D' and B' are ^uneration, dissipation and lateral boundary transports 
of perturbation kinetic energy, respectively. The generation term, due to the 
interaction of the zonal flow, the perturbation flow and the topography, is 


G' ■ - U £o f v* h ffl dx dy 
H J 

t>m + 0 


(3) 


* A 6ponge layer was introduced to dissipate waves with negative group 
velocity generated during the transient stage. Like the rigid top boundary 
condition of atmospheric models, the upstream boundary condition of unperturbed 
flow is well posed, but produces spurious reflection of transient long waves 
which propagate upstream. 
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Figure /. Open channel on a 0-plane geometry. The small sguad domain 
represents the region where pulsating lows are generated. 
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and is proportional to the Intensity of the northerly flow (v'<0) over the 
orography h m . This Indicates the existence of a positive feedback mechanism. 
Perturbation kinetic energy is realized, as Rossby lee waves , downstream of the 
mountain, whereas the upstream flow is largely undisturbed. Strong generation 
of perturbation kinetic energy will intensify the low in the lee of the mountain. 
This, in turn, Increases the zonal pressure gradient, increasing at the same 
time the northerly flow and the generation of perturbation energy. The positive 
feedback mechanism allows a simple interpretation of the multiple flow equi- 
libria solutions found by Charney and DeVore (1979) for nonlinear flow over 
topography. Under initial conditions corresponding to high index, the genera- 
tion of perturbation kinetic energy is inefficient, and the flow will remain 
in the "attractor basin" of the high index solution. If the Initial flow 
has a large amplitude perturbation flow (low index) with a favorable phase, 
there will be strong generation, and the flow will remain low index. If the 
phase is unfavorable perturbation kinetic energy will decay and the flow will 
evolve towards high index. 

The analysis of the generation term can be generalized to any type of stationary 
forcing (orography, heating or even locally enhanced /arocllnic instability). 

If we represent the effect of localized forcing by a term ^n(x,y) in the 
vorticity equation (1), the generation term becomes T 

G' - 1/t/J V^VF/dx dy (4) 

*0 + 0 


Again, the generation of perturbation energy depends on the phase relationship 
between the perturbation flow and the streamlines associated with the vorticity 
source. In the case of surface heating, there are also possible feedback 
mechanisms, since the intensity of the heat fluxes depends on the character- 
istics of the flow itself (speed, temperature, moisture content, etc.). 

In most of our experiments, we studied the interaction between stationary Rossby 
waves induced by orography (Figure 3a) and the flow induced by the periodic 
generation of pulses or "lows" upstream of the mountain. Every 3 days we 
increased the Instantaneous value of the potential vorticity by an amount 
equal to 30 percent of the Coriolis parameter. This was done in a square 
domain of 1500 km side indicated by darker shading in Figure 1. These "lows" 
can be interpreted as a crude simulation of the effect of periodic cyclogenesis 
in a region of enhanced baroclinicity. 

Since the pulses are always of the same sign, the net effect of the generation 
of lows in the absence of a mountain is to produce a small amplitude time- 
averaged flow (Figure 3b) superimposed with an even smaller amplitude time 
varying solution. It is the time-averaged solution that can Interact resonantly 
with the orography, so that our results are approximately independent of the 
frequency of generation of the pulses. 

When the pulses are generated upstream of the mountain, the results depend 
crucially on the position of the region of generation. If the pulses are 
generated at a region 6000 km upstream of the mountain, their induced flow is 
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Figure 3, Stream-function Field obtained with U ■ /Om sec~ } Ly - 5000 km, 
V = 20 days , maximum mountain height of 2000m, d) No pulses, 
b) Pulses (lows) every 3 days^ ho orography (time - averaged 
solution ); c and d ) Low s generated every 3 days From the. 
indicated sguare domain ; e) ,4.5 in 3b, but with maximum 
mountain height of 3000m , 
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unfavorable for the generation of perturbation kinetic energy over the 
orography the resulting flow (Figure 3c) is not very different from the 
unperturbed flow (Figure 3a). However, when the lows are generated at a 
region 4000 km upstream of the mountain, there is strong blocking development 
downstream of the orography (Figure 3d). Somewhat similar blocking patterns 
can be obtained by increasing the orographic forcing itself, from a maximum 
height of 2000 m to 3000 m (Figure 3e) in the absence of pulses. 

The orographic generation of perturbation kinetic energy corresponding to the 
three cases of Figure 3a, c and d is plotted in Figure 4. Note the very 
rapid response in the enhancement of the generation of KE' after the signal 
from the pulses first reach the orography. As indicated before, the generation 
is quite insensitive to the frequency of generation of individual pulses 
(except for its influence in the amplitude of the time-averaged flow). The 
observed overshooting in the generation is typical of the numerical experiments. 
The positive feedback mechanism increases the generation of KE' by deepening 
the low downstream of the mountain until the generation is compensated by 
advection and dissipation of KE' . This results, for smaller values of the 
dissipation constant (20 days in Figure 4), in a vacillating flow, even in 
the case of no pulses. Higher values of the friction ( T ■ 10 days) produce 
a steady state solution. 

The realistic blocking development presented in Figure 3d is a nonlinear 
resonance phenomenon. We repeated experiments 3a, 3c and 3d with a linear 
model (Figures 5a, 5b, and 5c). Even though there is, as expected, resonant 
enhancement, the pattern in Figure 5c does not resemble blocking. The effect 
of a rectified current downstream, referred to in Figure 2 and clearly ob- 
servable as a tilt in Figure 3 is absent in Figure 5. 

In Figure 6, we have increased the width of the channel from 5000 km, used 
in Figure 3 to 8000 km, so that stationary Rossby waves with meridional 
wavenumber 2 and 3 can also be excited. The results are similar to those 
previously presented, except that we observe a strong excitation of the 
meridional wavenumber 2 component downstream of the region of pulse genera- 
tion but immediately upstream of the mountain . Other integrations show an 
even stronger wavenumber 2 response, resembling the observed blocking upstream 
of the Rockies with low pressure in the coast of California and an intense 
ridge in northern latitudes. 

The simple model presented here may be relevant to blocking in the real atmos- 
sphere. In particular, persistent blocking in the Atlantic seems associated 
with a favorable circulation over the Rockies, with enhanced northerly flow 
over orography and perhaps enhanced heating over the Gulf Stream. This can 
be seen to be the case in the persistent anomalous circulation observed in 
December 1976 (Taubensee, 1977) and January 1977 (Wagner, 1977). Indeed, 
enhanced orographic forcing seems to have been associated with a favorable 
distribution of anomalous sea surface temperatures in the Pacific. Following 
a suggestion of Namias (1978), Shukla and Bangaru (1979) have recently per- 
formed a sensitivity experiment with a GLAS GCM in which they introduced 
observed sea surface temperature anomalies in the Pacific multiplied by a 
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Figure, S, As in Figure s 3a, 3c and 3d , bu1~ wiih 0 linear mode./. 


1 


f 





285 


J 



factor of 9/5 (Figure 7a). The SST anomalies enhanced the generation term 
over the Rockies (Figure 7b) and, as in the simple barotropic model, produced 
a maximum response downstream of the regions of orographic and thermal forcing. 

We are presently conducting a systematic study of observed Atlantic blocking 
situations, and the relationship between the generation of perturbation energy 
over the Rockies and the duration of blocking. 
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Figure 7 a. Observed Pacific Sea Surface Temperature Anomalies (°F) durino 
the Winter of 1976- 1977. * 
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